r7.01.32 Loi de comportement de plaques en béton armé GLRC_DM

Résumé :

Cette documentation présente la formulation théorique et l’intégration numérique de la loi de comportement GLRC_DM, utilisable avec la modélisation DKTG. Elle fait partie des modèles dits «globaux» utilisés pour des structures minces (poutres, plaques et coques). Les phénomènes non-linéaires, tels que la plasticité ou l’endommagement, sont directement en relations avec les déformations généralisées (extension, courbure, distorsion) et les contraintes généralisées (efforts de membrane, de flexion et tranchants). Ainsi, cette loi de comportement s’applique avec un élément fini de plaque ou de coque. Cela permet d’économiser, par rapport à une approche multi-couches, du temps CPU ainsi que de la mémoire. L’avantage par rapport aux coques multi-couches est encore plus important, lorsqu’un des constituants de la plaque se comporte de manière quasi-fragile (béton, par exemple), puisque le modèle global permet d’éviter les problèmes de localisation.

La loi de comportement GLRC_DM modélise l’endommagement sous sollicitation membranaire et sous sollicitation de flexion de plaques en béton armé, à l’aide de paramètres «homogénéisés». Ce modèle de comportement représente donc une évolution par rapport au modèle GLRC_DAMAGE qui ne traite l’endommagement qu’en sollicitation de flexion. La structure du modèle d’endommagement de GLRC_DM ressemble à celle de GLRC_DAMAGE et sont toutes les deux inspirées de ENDO_ISOT_BETON. Contrairement à GLRC_DAMAGE, GLRC_DM ne permet pas de modéliser une éventuelle plastification, ce qui le rend moins adapté pour simuler des chargements extrêmes. De plus, le comportement GLRC_DM est isotrope avant endommagement: on néglige l’orthotropie apportée par les nappes d’armatures en acier.

On peut envisager d’utiliser le comportement GLRC_DM «seul» sur une modélisation de plaques pour représenter le béton armé, ou bien uniquement pour représenter le béton seulement en l’associant alors avec des modélisations de grilles d’armatures d’acier, ce qui permet de représenter l’orthotropie et éventuellement la non symétrie des nappes d’armatures. Ce dernier choix simplifie aussi le paramétrage de ce comportement.

Formulation du modèle de comportement

L’emploi de la théorie des plaques et coques minces permet de décrire efficacement le comportement mécanique des structures en béton armé, qui sont généralement élancées.

Dans une première étape de la construction du modèle, on suppose l’existence d’un milieu homogénéisé avec le même comportement mécanique que la structure en béton armé, à laquelle on s’intéresse. Pour simplifier, on fait l’hypothèse que ce milieu est isotrope et que l’élément de structure étudié est symétrique par rapport à son feuillet moyen. Ces hypothèses ne sont pas indispensables pour la formulation, mais ont été faites pour simplifier la démarche. De plus, on estime que leur impact sur le comportement est moindre par rapport à la fissuration, qui est au centre d’intérêt du modèle.

On doit noter que l’utilisation de ce modèle est associée à celle d’un élément de plaque. Il n’est utilisable que dans le cadre d’éléments finis DKT (modélisation supportée : DKTG), correspondant à la théorie de Love-Kirchhoff , où l’on néglige toute distorsion transverse dans l’épaisseur de la plaque.

Énergie libre

Pour un milieu continu homogène isotrope au comportement élastique linéaire on peut écrire la densité volumique de l’énergie libre comme:

\({\Phi}_{e}(\varepsilon )=\frac{\lambda}{2}\mathrm{tr}{(\varepsilon )}^{2}+\mu \sum_{i=1}^{3}{{\tilde{\varepsilon}}_{i}}^{2}\)

où \(\lambda ,\mu\) sont les coefficients de Lamé, \(\varepsilon\) le tenseur de déformation et \({\tilde{\varepsilon}}_{i}\) ses valeurs propres. Comme pour le modèle ENDO_ISOT_BETON on introduit l’endommagement par une fonction \(\xi (.,{d}_{i})\) , \({d}_{i}\) étant une variable d’endommagement. Donc, pour un milieu endommageable, l’énergie libre s’écrit:

(3675)\[{\Phi}_{\mathrm{ed}}(\varepsilon ,{d}_{j})=\frac{\lambda}{2}\mathrm{tr}{(\varepsilon )}^{2}\xi (\mathrm{tr}(\varepsilon ),{d}_{j})+\mu \sum_{i=1}^{3}{{\tilde{\varepsilon}}_{i}}^{2}\xi ({\tilde{\varepsilon}}_{i},{d}_{j})\]

avec la fonction \((x,d)\in {ℝ}^{2}\to \xi (x,d)\le 1\) vérifiant a priori \(\frac{\partial \xi }{\partial d}\le 0\) pour représenter la perte de raideur liée à l’endommagement, et \(\frac{\partial \xi }{\partial x}=0 \quad \mathrm{pour} \quad x\in \left]-\infty ,0\right[\cup \left]0,+\infty \right[\) , le saut en \(0\) permettant de représenter la discontinuité de comportement entre traction et compression.

L’équation (4894) est valable pour un milieu continu et on va l’appliquer à une plaque de Love-Kirchhoff \(\omega \times \left]\frac{-h}{2};\frac{h}{2}\right[\) , d’épaisseur \(h\) (on note \(z={x}_{3}\) ), où on fait des hypothèses cinématiques de Hencky-Mindlin (voir [bib9]):

\(\left(\begin{array}{}{U}_{1}({x}_{1,}{x}_{2,}z)\\ {U}_{2}({x}_{1,}{x}_{2,}z)\\ {U}_{z}({x}_{1,}{x}_{2,}z)\end{array} \right)=\underset{\begin{array}{}\mathrm{cinématique}\mathrm{de}\mathrm{plaque}\\ {u}^{s}\in {V}_{s}\end{array}}{\underset{\underbrace{}}{\left(\begin{array}{}{u}_{1}({x}_{1,}{x}_{2})\\ {u}_{2}({x}_{1,}{x}_{2})\\ {u}_{z}({x}_{1,}{x}_{2})\end{array} \right)+z\left(\begin{array}{}{\theta}_{2}({x}_{1,}{x}_{2})\\ -{\theta}_{1}({x}_{1,}{x}_{2})\\ 0\end{array} \right)}}+\underset{\begin{array}{}\mathrm{déplacement}\mathrm{complémentaire}\\ {u}^{c}\in {V}_{c}\end{array}}{\underset{\underbrace{}}{\left(\begin{array}{}{u}_{1}^{c}({x}_{1,}{x}_{2,}z)\\ {u}_{2}^{c}({x}_{1,}{x}_{2,}z)\\ {u}_{z}^{c}({x}_{1,}{x}_{2,}z)\end{array} \right)}}\)

où \(U={({U}_{1}{U}_{2}{U}_{z})}^{T}\) est le champ de déplacement en 3D, \(u={({u}_{1}{u}_{2}{u}_{z})}^{T}\) le déplacement du feuillet moyen et \({\theta}_{x}\) , \({\theta}_{y}\) ses rotations. Ainsi, le tenseur de déformation, défini comme:

\({\varepsilon}_{ij}=\frac{1}{2}\left(\frac{\partial {U}_{i}}{\partial {x}_{j}}+\frac{\partial {U}_{j}}{\partial {x}_{i}}\right),i,j=\mathrm{1..3}\)

s’écrit aussi comme:

(3676)\[\begin{split}\begin{array}{}{\varepsilon}_{11}=\underset{{\varepsilon}_{11}^{s}}{\underset{\underbrace{}}{{\epsilon}_{11}+z{\kappa}_{11}}}+{u}_{1,1}^{c}\\ {\varepsilon}_{22}=\underset{{\varepsilon}_{22}^{s}}{\underset{\underbrace{}}{{\epsilon}_{22}+z{\kappa}_{22}}}+{u}_{2,2}^{c}\\ {\varepsilon}_{12}=\underset{{\varepsilon}_{12}^{s}}{\underset{\underbrace{}}{{\epsilon}_{12}+z{\kappa}_{12}}}+\frac{1}{2}({u}_{2,1}^{c}+{u}_{1,2}^{c})\\ {\varepsilon}_{\mathrm{1z}}={\varepsilon}_{\mathrm{1z}}^{c}=\frac{1}{2}{u}_{3,1}^{c}\\ {\varepsilon}_{\mathrm{2z}}={\varepsilon}_{\mathrm{2z}}^{c}=\frac{1}{2}{u}_{3,2}^{c}\\ {\varepsilon}_{zz}={\varepsilon}_{zz}^{c}={u}_{3,z}^{c}\end{array}\end{split}\]

où \(\epsilon\) est le tenseur de l’extension membranaire, défini dans le plan:

\({\epsilon}_{ij}=\frac{1}{2}\left(\frac{\partial {u}_{i}}{\partial {x}_{j}}+\frac{\partial {u}_{j}}{\partial {x}_{i}}\right),i,j=\mathrm{1..2}\)

et \(\kappa\) le tenseur de variation de courbure, défini dans le plan:

\({\kappa}_{11}=\frac{\partial {\theta}_{2}}{\partial {x}_{1}}\) , \({\kappa}_{22}=\frac{-\partial {\theta}_{1}}{\partial {x}_{2}}\) , \({\kappa}_{12}=\frac{1}{2}\left(\frac{\partial {\theta}_{2}}{\partial {x}_{2}}-\frac{\partial {\theta}_{1}}{\partial {x}_{1}}\right)\)

relations auxquelles se rajoute l’hypothèse de contraintes planes \({\sigma}_{zz}=0\) , \({\sigma}_{\mathrm{1z}}=0\) , \({\sigma}_{\mathrm{2z}}=0\) qui déterminera le champ de déplacement complémentaire \({u}^{c}\in {V}_{c}\) . Dans la théorie utilisée ici, on n’introduit que deux composantes de rotation \({\theta}_{x}\) et \({\theta}_{y}\) , ce qui implique que le tenseur de variation de courbure est 2D et n’a que 3 composantes indépendantes.

En introduisant ces hypothèses cinématiques, cf. (4895), dans l’expression de l’énergie libre, (4894), on peut déterminer les valeurs propres de la déformation \({\tilde{\varepsilon}}_{i}=\tilde{{(\epsilon +z\kappa )}_{i}}\) . Ces valeurs propres étant, en général, des fonctions non polynomiales de la coordonnée \(z\) , l’intégrale \(\int{\phi }_{\mathrm{ed}}\mathrm{dz}\) n’est pas calculable analytiquement. Cette formulation ne serait donc pas adaptée pour l’application aux éléments de structure.

Au lieu de la formulation (4894), on utilisera plutôt une formulation de l’énergie libre définie directement en déformations généralisées, \(\epsilon\) et \(\kappa\) :

(3677)\[\begin{split}\begin{array}{ccc}{\Phi}_{\mathrm{ed}}(\epsilon ,\kappa ,{\varepsilon}_{zz},{d}_{j})& =& \frac{\lambda}{2}{(\mathrm{tr}(\epsilon )+{\varepsilon}_{zz})}^{2}.{\xi}_{m}(\mathrm{tr}(\epsilon ),{d}_{j})+\mu (\sum_{i=1}^{2}{{\tilde{\epsilon}}_{i}}^{2}.{\xi}_{m}({\tilde{\epsilon}}_{i},{d}_{j})+{\varepsilon}_{zz}^{2})\\ & +& \frac{\lambda}{2}{z}^{2}\mathrm{tr}{(\kappa )}^{2}.{\xi}_{f}(\mathrm{tr}(\kappa ),{d}_{j})+\mu {z}^{2}\sum_{i=1}^{2}{{\tilde{\kappa}}_{i}}^{2}.{\xi}_{f}({\tilde{\kappa}}_{i},{d}_{j})+z(.)\end{array}\end{split}\]

où \(z(.)\) contient tous les termes couplant \(\epsilon\) et \(\kappa\) , qui disparaissent après l’intégration sur \(z\) , si l’on fait l’hypothèse que la plaque/poutre est symétrique par rapport au feuillet moyen. On obtient ainsi la densité surfacique de l’énergie libre:

(3678)\[\begin{split}\begin{array}{ccc}{\Phi}_{\mathrm{ed}}^{S}(\epsilon ,\kappa ,{\varepsilon}_{zz},{d}_{j})& =& {\int}_{-\frac{h}{2}}^{\frac{h}{2}}{\Phi}_{\mathrm{ed}}(\epsilon ,\kappa ,{\varepsilon}_{zz},{d}_{j})\mathrm{dz}\\ & =& \frac{{\lambda}_{m}}{2}{(\mathrm{tr}(\epsilon )+{\varepsilon}_{zz})}^{2}.{\xi}_{m}(\mathrm{tr}(\epsilon ),{d}_{j})+{\mu}_{m}(\sum_{i=1}^{2}{{\tilde{\epsilon}}_{i}}^{2}.{\xi}_{m}({\tilde{\epsilon}}_{i},{d}_{j})+{\varepsilon}_{zz}^{2})\\ & +& \frac{{\lambda}_{f}}{2}\mathrm{tr}{(\kappa )}^{2}.{\xi}_{f}(\mathrm{tr}(\kappa ),{d}_{j})+{\mu}_{f}\sum_{i=1}^{2}{{\tilde{\kappa}}_{i}}^{2}.{\xi}_{f}({\tilde{\kappa}}_{i},{d}_{j})\end{array}\end{split}\]

où \(h\) est l’épaisseur de la plaque et \({\lambda}_{m},{\lambda}_{f},{\mu}_{m},{\mu}_{f}\) sont définis à la section r7.01.32-parametres-elastique-endommageable . La densité surfacique de l’énergie libre est strictement convexe pour \({d}_{j}=0\) (ayant correctement choisi les coefficients élastiques).

Dans (4896) et (4897) nous avons également fait l’hypothèse que l’endommagement n’est pas affecté par l’extension en \(z\) , ce qui se traduit par l’absence de \({\varepsilon}_{zz}\) des arguments de l’indicateur d’endommagement \({\xi}_{m}\) . Cela se justifie par l’objectif de ce modèle qui est de représenter la fissuration perpendiculaire par rapport au feuillet moyen, qui se déclenche soit en sollicitation membranaire soit en flexion, mais jamais en extension par \({\varepsilon}_{zz}\) . De plus, du point de vue de résolution numérique abordée plus bas, cette hypothèse facilite le calcul local des variables d’endommagement \({d}_{j}\) et la satisfaction de l’hypothèse de contrainte plane \({\sigma}_{zz}=0\) .

Remarque

Dans (4897), on observe que la déformation* \({\varepsilon}_{zz}\) n’est introduite explicitement que sur le terme d’énergie membranaire. L’effet de contraintes planes sur le comportement membranaire sera donc affecté par l’endommagement. Ceci n’est pas le cas en flexion. On pourrait imaginer que la condition de contraintes planes soit écrite en intégrant un couplage flexion-membrane, lui-même dépendant de l’endommagement, mais cela n’a pas été chois ici.

On verra section r7.01.32-parametres-elastique-endommageable comment déterminer les paramètres du modèle à partir de ses réponses sur des cas simples.

Quant à la variable d’endommagement, elle est décrite par deux composantes, l’une des deux représentant globalement l’endommagement du côté de la face supérieure de la plaque et l’autre représentant l’endommagement du côté de la face inférieure de la plaque:

\(d(z)=\left\{ \begin{array}{ccc}{d}_{1}& \text{si}& z\ge 0\\ {d}_{2}& \text{si}& z<0\end{array}\right.\)

Il reste à définir les fonctions caractéristiques de l’endommagement \({\xi}_{m}(.,{d}_{j})\) et \({\xi}_{f}(.,{d}_{j})\) , pour que la formulation, (4897), soit complète:

\({\xi}_{m}(x,{d}_{1,}{d}_{2})=\frac{1}{2}((\frac{1+{\gamma}_{\mathrm{mt}}{d}_{1}}{1+{d}_{1}}+\frac{1+{\gamma}_{\mathrm{mt}}{d}_{2}}{1+{d}_{2}})H(x)+(\frac{{\alpha}_{c}+{\gamma}_{\mathrm{mc}}{d}_{1}}{{\alpha}_{c}+{d}_{1}}+\frac{{\alpha}_{c}+{\gamma}_{\mathrm{mc}}{d}_{2}}{{\alpha}_{c}+{d}_{2}})H(-x))\in [0,1]\)

et

\({\xi}_{f}(x,{d}_{1,}{d}_{2})=\frac{\alpha +{\gamma}_{f}{d}_{1}}{\alpha +{d}_{1}}H(x)+\frac{\alpha +{\gamma}_{f}{d}_{2}}{\alpha +{d}_{2}}H(-x)\in [0,1]\)

où \(H(.)\) est la fonction de Heaviside.

Les fonctions \((x,d)\in {ℝ}^{2}\to \xi (x,d)\le 1\) , qui s’expriment de la même manière que celle choisie dans la loi ENDO_ISOT_BETON, cf. [r7.01.04], en opérant un changement de variable sur \({d}_{i}\) , offrent l’avantage de donner une pente constante dans les phases d’endommagement. Elles sont décroissantes \(\frac{\partial \xi }{\partial d}\le 0\) et convexes pour représenter la perte de raideur liée à l’endommagement, et \(\frac{\partial \xi }{\partial x}=0\mathrm{pour}x\in \left]-\infty ,0\right[\cup \left]0,+\infty \right[\) , le saut en \(0\) permettant de représenter le changement de comportement entre traction et compression, cependant sans introduire de discontinuité. Les variables d’endommagement \({d}_{i}\) croissent jusqu’à \(+\infty\) .

Les fonctions caractéristiques de l’endommagement \({\xi}_{m}(.,{d}_{j})\) et \({\xi}_{f}(.,{d}_{j})\) varient de \(1\) à respectivement \({\gamma}_{\mathrm{mt}}\) ou \({\gamma}_{\mathrm{mc}}\) , et \({\gamma}_{f}\) , pour des \({d}_{j}\to +\infty\) . Ce modèle décrit donc un endommagement partiel mais pas total du matériau.

Les paramètres d’endommagement \({\gamma}_{\mathrm{mt}}\) pour la traction en membrane, \({\gamma}_{\mathrm{mc}}\) pour la compression en membrane et \({\gamma}_{f}\) pour la flexion, peuvent avoir des valeurs dans \([0,1]\) , de manière à ce que le modèle ne soit pas adoucissant, ce qui entraînerait des difficultés de dépendance à la discrétisation spatiale et de convergence. On choisira \({\gamma}_{.}\approx 0\) lorsque le phénomène correspondant aura plus d’impact sur l’endommagement et \({\gamma}_{.}\approx 1\) lorsque celui-ci sera négligeable. Ainsi, pour le béton armé, on s’attend à \({\gamma}_{\mathrm{mc}}\approx 1\) et \({\gamma}_{\mathrm{mt}}\approx 0\) . Quant au paramètre, \(\alpha\) il permet d’ajuster la contribution de la flexion au seuil d’endommagement. Le paramètre \({\alpha}_{c}\) permet de moduler l’évolution de l’endommagement en compression.

Contraintes généralisées

Selon la procédure habituelle on définit les contraintes généralisées (efforts normaux et moments) par les dérivées de la densité d’énergie libre par rapport aux déformations généralisées, \(\epsilon\) et \(\kappa\) :

\(N=\frac{\partial {\Phi}_{\mathrm{ed}}^{S}}{\partial \epsilon }\) ; \(M=\frac{\partial {\Phi}_{\mathrm{ed}}^{S}}{\partial \kappa }\)

Dans notre application, les contraintes généralisées sont calculées dans le repère des vecteurs propres des déformations généralisées et elles s’écrivent comme (voir le calcul en annexe):

(3679)\[\begin{split}\begin{array}{ccccc}{\tilde{N}}_{i}& =& \frac{\partial {\Phi}_{\mathrm{ed}}^{S}}{\partial {\tilde{\epsilon}}_{i}}=& {\lambda}_{m}(\mathrm{tr}(\epsilon )+{\varepsilon}_{zz}).{\xi}_{m}(\mathrm{tr}(\epsilon ),{d}_{j})+2{\mu}_{m}{\tilde{\epsilon}}_{i}.{\xi}_{m}({\tilde{\epsilon}}_{i},{d}_{j}),& i=\mathrm{1..2}\\ {\tilde{M}}_{i}& =& \frac{\partial {\Phi}_{\mathrm{ed}}^{S}}{\partial {\tilde{\kappa}}_{i}}=& {\lambda}_{f}\mathrm{tr}(\kappa ).{\xi}_{f}(\mathrm{tr}(\kappa ),{d}_{j})+2{\mu}_{f}{\tilde{\kappa}}_{i}.{\xi}_{f}({\tilde{\kappa}}_{i},{d}_{j}),& i=1,2\end{array}\end{split}\]

Il est facile de vérifier (voir annexe) que les contraintes généralisées \(N\) et \(M\) , définies à partir des expressions (4898), partagent respectivement les mêmes vecteurs propres que les déformations généralisées \(\epsilon\) et \(\kappa\) . Si ces vecteurs propres sont désignés respectivement comme \({Q}_{m}\) et \({Q}_{f}\) , on peut écrire:

(3680)\[N={Q}_{m}\tilde{N}{Q}_{m}^{T} \quad \text{et} \quad M={Q}_{f}\tilde{M}{Q}_{f}^{T}\]

où \(\tilde{N}\) et \(\tilde{M}\) sont les matrices diagonales composées des valeurs propres définies dans (4898). Il est important de noter que les vecteurs propres pour la partie de membrane \({Q}_{m}\) et de flexion et \({Q}_{f}\) sont complètement indépendants.

De même, on définit les contraintes de pincement comme la variable duale de \({\varepsilon}_{zz}\) :

(3681)\[{\sigma}_{zz}=\frac{\partial {\Phi}_{\mathrm{ed}}^{S}}{\partial {\varepsilon}_{zz}}={\lambda}_{m}{\xi}_{m}(\mathrm{tr}(\epsilon ),{d}_{1,}{d}_{2}).(\mathrm{tr}(\epsilon )+{\varepsilon}_{zz})+2{\mu}_{m}{\varepsilon}_{zz}\]

à laquelle on imposera la condition de contraintes planes: \({\sigma}_{zz}=0\) .

Seuils et évolution de l’endommagement

Pour pouvoir définir un seuil d’endommagement dans le cadre de l’hypothèse d’un matériau standard généralisé (voir [bib1], [bib7]), on introduit les forces thermodynamiques associées aux variables \({d}_{1}\) et \({d}_{2}\) :

(3682)\[{Y}_{j}=-\frac{\partial {\Phi}_{\mathrm{ed}}^{S}}{\partial {d}_{j}}={Y}_{j}^{m}+{Y}_{j}^{f}\]

où

(3683)\[ \begin{align}\begin{aligned}{Y}_{j}^{m}(\epsilon ,{d}_{j},{\varepsilon}_{zz})=\frac{1}{{(1+{d}_{j})}^{2}}(\frac{{\lambda}_{m}}{4}{(\mathrm{tr}(\epsilon )+{\varepsilon}_{zz})}^{2}.{G}_{m}(\mathrm{tr}(\epsilon ),{d}_{j})+\frac{{\mu}_{m}}{2}(\sum_{i=1}^{2}{\tilde{\epsilon}}_{i}^{2}.{G}_{m}({\tilde{\epsilon}}_{i},{d}_{j})))\\\text{en notant} \quad {G}_{m}(x,{d}_{j})=(1-{\gamma}_{\mathrm{mt}})H(x)+(\frac{{\alpha}_{c}(1-{\gamma}_{\mathrm{mc}}){(1+{d}_{j})}^{2}}{{({\alpha}_{c}+{d}_{j})}^{2}})H(-x)\in [0,1]\\\text{de telle sorte que :} \quad \frac{\partial {\xi}_{m}(x,{d}_{1,}{d}_{2})}{\partial {d}_{j}}=-\frac{{G}_{m}(x,{d}_{j})}{2{(1+{d}_{j})}^{2}}\\\text{et} \quad {Y}_{j}^{f}(\kappa ,{d}_{j})=\frac{\alpha}{{(\alpha +{d}_{j})}^{2}}{Y}_{j}^{f,0}(\kappa )\\\text{avec} \quad {Y}_{j}^{f,0}(\kappa )=(1-{\gamma}_{f})(\frac{{\lambda}_{f}}{2}\mathrm{tr}{(\kappa )}^{2}H({(-1)}^{(j+1)}.\mathrm{tr}(\kappa ))+{\mu}_{f}\sum_{i=1}^{2}{\tilde{\kappa}}_{i}^{2}H({(-1)}^{(j+1)}.{\tilde{\kappa}}_{i}))\end{aligned}\end{align} \]

Les seuils d’endommagement sont définis par:

(3684)\[{f}_{{d}_{j}}={Y}_{j}(\epsilon ,\kappa ,{d}_{j},{\varepsilon}_{zz})-{k}_{{0}_{j}}\le 0\]

où \({k}_{{0}_{j}}\) sont des constantes de seuil. Ces seuils définissent le domaine convexe de réversibilité dans l’espace des \((\epsilon ,\kappa )\) .

En principe les constantes de seuil \({k}_{{0}_{1}}\) et \({k}_{{0}_{2}}\) pourraient être différentes, mais selon l’hypothèse qu’on fait sur la symétrie de la plaque par rapport au feuillet moyen, les deux valeurs sont les mêmes: \({k}_{{0}_{1}}={k}_{{0}_{2}}={k}_{0}\) .

On voit dans (3683) que le paramètre \(\alpha\) régule la contribution de la flexion au seuil d’endommagement initial, puisque:

\({f}_{{d}_{j}}({d}_{j}=0)={Y}_{j}^{m}(\epsilon ,{d}_{j}=0)+\frac{1}{\alpha}{Y}_{j}^{f,0}(\kappa )-{k}_{0}\)

La loi d’évolution des variables d’endommagement \({d}_{1}\) et \({d}_{2}\) est définie par la règle de normalité aux seuils (3684), pour lesquels on peut définir le pseudo-potentiel de dissipation \(D(\delta )\) :

(3685)\[\begin{split}\begin{array}{}{\dot{d}}_{j}=\eta \frac{\partial {f}_{{d}_{j}}}{\partial {Y}_{j}},\mathrm{avec}\eta \ge 0\\ \iff {Y}_{j}\in \partial D({\dot{d}}_{j})\iff D({\dot{d}}_{j})-D(\delta )\ge {Y}_{\mathrm{j.}}({\dot{d}}_{j}-\delta ),\forall \delta \ge 0\end{array}\end{split}\]

Les valeurs d’endommagement \({d}_{1}\) et \({d}_{2}\) sont parfaitement déterminées par les conditions suivantes:

(3686)\[\textrm{si} {f}_{{d}_{j}}<0 \textrm{si} {f}_{{d}_{j}}=0\]

L’évolution des variables d’endommagement est donc obtenue à l’aide de la condition de cohérence, les fonctions \(\xi (x,d)\) étant convexes, les modules d’écrouissage \(-{f}_{,{d}_{j}}\) sont positifs (les coefficients vérifiant \({\gamma}_{.}\in [0,1]\) ):

(3687)\[{\dot{d}}_{j}=-\frac{{[{f}_{,Y}.{Y}_{,\epsilon }.\dot{\epsilon}]}_{+}}{{f}_{,{d}_{j}}}\]

On constate qu’en charge membranaire pure uniaxiale endommageante (\(\dot{{d}_{1}}=\dot{{d}_{2}}\ge 0\) ), la force thermodynamique s’exprime: \({Y}^{m}(d)=-{E}_{\mathrm{eq}}^{m}h\frac{{\varepsilon}^{2}.{\xi}_{m,d}(d)}{2}\) , \({E}_{\mathrm{eq}}^{m}h\) étant la raideur élastique membranaire uniaxiale. La condition de cohérence s’écrit alors:

(3688)\[\dot{{f}_{d}}={E}_{\mathrm{eq}}^{m}h\left(-\epsilon \dot{\epsilon}{\xi}_{m,d}-\frac{{\epsilon}^{2}\dot{d}{\xi}_{m,\mathrm{dd}}}{2}\right)=0\]

D’où: \(\dot{{d}_{m}}=-\frac{2\dot{\epsilon}{\xi}_{m,d}}{\epsilon {\xi}_{m,\mathrm{dd}}}\) et donc la loi de charge membranaire pure uniaxiale endommageante est:

(3689)\[\dot{N}={E}_{\mathrm{eq}}^{m}h({\xi}_{m,d}\dot{d}\epsilon +{\xi}_{m}\dot{\epsilon})={E}_{\mathrm{eq}}^{m}h\dot{\epsilon}\left({\xi}_{m}-2\frac{{\xi}_{m,d}^{2}}{{\xi}_{m,\mathrm{dd}}} \right)={\gamma}_{m}{E}_{\mathrm{eq}}^{m}h\dot{\epsilon}\]

ce qui permet d’interpréter le rôle du paramètre \({\gamma}_{m}\) : la pente étant constante, ce qui est la justification de la forme algébrique de la fonction \({\xi}_{m}\) .

Intégration numérique

Contrairement à la majorité des lois de comportement anélastiques non linéaires, celle présentée ici ne nécessite pas la discrétisation des équations d’évolution des variables internes. On adopte une méthode de discrétisation directe implicite en temps.

Les variables \({d}_{j}\) peuvent être calculées directement à partir de la condition de cohérence, (3685). La seule fois où on fait référence à la «vitesse d’endommagement» c’est pour vérifier qu’elle est positive \({\dot{d}}_{j}\ge 0\) . Pour un calcul incrémental cette condition se traduit par \({d}_{j}^{n}\ge {d}_{j}^{n-1}\) au pas de temps \(n>1\) , donc sans référence à un schéma particulier d’intégration temporelle.

Plaçons à un instant donné \({t}_{n}\) du trajet de chargement. On effectue d’abord une étape de prédiction élastique (tenseur d’élasticité évalué avec les variables d’endommagement \({d}_{j}^{n-1}\) figées à l’étape précédente), d’où \(({\epsilon}^{n},{\kappa}^{n})\) . On calcule ensuite un nouveau \({\varepsilon}_{zz}^{\mathrm{n0}}\) :

\({\sigma}_{zz}^{n}=0\Rightarrow {\varepsilon}_{zz}^{\mathrm{n0}}({\epsilon}^{n},{d}_{1}^{n-1},{d}_{2}^{n-1})=-\frac{{\lambda}_{m}{\xi}_{m}(\mathrm{tr}({\epsilon}^{n}),{d}_{1}^{n-1},{d}_{2}^{n-1})}{2{\mu}_{m}+{\lambda}_{m}{\xi}_{m}(\mathrm{tr}({\epsilon}^{n}),{d}_{1}^{n-1},{d}_{2}^{n-1})}\mathrm{tr}({\epsilon}^{n})\)

voir (3681).

On calcule ensuite \({f}_{{d}_{j}}^{\mathrm{n0}}({\epsilon}^{n},{\kappa}^{n},{d}_{j}^{n-1},{\varepsilon}_{zz}^{\mathrm{n0}})\) pour vérifier les seuils d’endommagement. Si \({f}_{{d}_{j}}^{\mathrm{n0}}\le 0\) , l’endommagement n’évolue pas: \({d}_{j}^{n}={d}_{j}^{n-1}\) et \({\varepsilon}_{zz}^{n}={\varepsilon}_{zz}^{\mathrm{n0}}\) , et les contraintes généralisées se calculent selon (4898).

Lorsque \({f}_{{d}_{j}}^{\mathrm{n0}}>0\) , l’endommagement peut évoluer et on doit résoudre les équations:

\({f}_{{d}_{j}}^{n}={Y}_{j}({\epsilon}^{n},{\kappa}^{n},{d}_{j}^{n},{\varepsilon}_{zz}^{n})-{k}_{0}=0\)

Ce qui correspond à la résolution des équations non linéaires, à \(\epsilon ,\kappa\) donnés:

(3690)\[\begin{split}\begin{array}{}{R}_{{d}_{j}}({d}_{j},\epsilon ,\kappa ,{\varepsilon}_{zz})=\frac{1}{{(1+{d}_{j})}^{2}}(\frac{{\lambda}_{m}}{4}{(\mathrm{tr}(\epsilon )+{\varepsilon}_{zz})}^{2}.{G}_{m}(\mathrm{tr}(\epsilon ),{d}_{j})+\frac{{\mu}_{m}}{2}(\sum_{i=1}^{2}{\tilde{\epsilon}}_{i}^{2}{G}_{m}({\tilde{\epsilon}}_{i},{d}_{j})))+\\ \frac{\alpha (1-{\gamma}_{f})}{{(\alpha +{d}_{j})}^{2}}(\frac{{\lambda}_{f}}{2}\mathrm{tr}{(\kappa )}^{2}H({(-1)}^{(j+1)}\mathrm{tr}(\kappa ))+{\mu}_{f}\sum_{i=1}^{2}{\tilde{\kappa}}_{i}H({(-1)}^{(j+1)}{\tilde{\kappa}}_{i}))-{k}_{0}=0\end{array}\end{split}\]

L’équation (3690) doit être résolue en prenant en compte également la condition de contrainte plane, qui permet d’exprimer \({\varepsilon}_{zz}(\epsilon ,{d}_{j})\) , cf. (3681):

(3691)\[{\sigma}_{zz}=0\Rightarrow {\varepsilon}_{zz}(\epsilon ,{d}_{1,}{d}_{2})=-\frac{{\lambda}_{m}{\xi}_{m}(\mathrm{tr}(\epsilon ),{d}_{1,}{d}_{2})}{2{\mu}_{m}+{\lambda}_{m}{\xi}_{m}(\mathrm{tr}(\epsilon ),{d}_{1,}{d}_{2})}\mathrm{tr}(\epsilon )\]

On résout les équations (3690) et (3691) par la méthode de Newton. On commence par une phase de prédiction, de type Euler explicite:

\({d}_{j}^{n(0)}={d}_{j}^{n-1}+{\left(\frac{{\mathrm{dR}}_{{d}_{j}}}{d({d}_{j})} \right)\mid }_{{d}_{j}^{n-1}}^{-1}.({k}_{0}-{Y}_{j}({\epsilon}^{n},{\kappa}^{n},{d}_{j}^{n-1},{\varepsilon}_{zz}^{\mathrm{n0}}))\)

Puis on traite la phase de correction, à l’itération \((m)>1\) :

1. \(\Delta {d}_{j}^{n(m)}={\left(\frac{{\mathrm{dR}}_{{d}_{j}}}{d({d}_{j})}\right)\mid }_{{d}_{j}^{n(m-1)}}^{-1}.({k}_{0}-{Y}_{j}({\epsilon}^{n},{\kappa}^{n},{d}_{j}^{n(m-1)},{\varepsilon}_{zz}^{n(m-1)}))\)

2. \({d}_{j}^{n(m)}={d}_{j}^{n(m-1)}+\Delta {d}_{j}^{n(m)}\)

3. \({\varepsilon}_{zz}^{n(m)}={\varepsilon}_{zz}({\epsilon}^{n},{d}_{1}^{n(m)},{d}_{2}^{n(m)})\)

Cette phase de correction s’achève lorsque le critère de convergence exprimé en terme énergétique est atteint:

\(\Delta {d}_{j}^{n(m)}.{R}_{{d}_{j}}^{n(m)}<{\eta}_{\mathrm{tolerance}}.{k}_{0}\)

L’opérateur tangent de ce système non linéaire est défini comme:

(3692)\[\begin{split}\begin{array}{ccc}\frac{{\mathrm{dR}}_{{d}_{j}}}{d({d}_{j})}& =& \frac{\partial {R}_{{d}_{j}}}{\partial {\varepsilon}_{zz}}\frac{\partial {\varepsilon}_{zz}}{\partial {d}_{j}}+\frac{\partial {R}_{{d}_{j}}}{\partial {d}_{j}}\\ & =& \frac{1}{{(1+{d}_{j})}^{2}}(\frac{{\lambda}_{m}}{2}(\mathrm{tr}(\epsilon )+{\varepsilon}_{zz}){G}_{m}(\mathrm{tr}(\epsilon ),{d}_{j})\frac{\partial {\varepsilon}_{zz}}{\partial {d}_{j}})\\ & +& \frac{1}{{(1+{d}_{j})}^{2}}(\frac{{\lambda}_{m}}{4}{(\mathrm{tr}(\epsilon )+{\varepsilon}_{zz})}^{2}\frac{\partial {G}_{m}(\mathrm{tr}(\epsilon ),{d}_{j})}{\partial {d}_{j}}{\varepsilon}_{zz}+\frac{{\mu}_{m}}{2}(\sum_{i=1}^{2}{\tilde{\epsilon}}_{i}^{2}\frac{\partial {G}_{m}({\tilde{\epsilon}}_{i},{d}_{j})}{\partial {d}_{j}}))\\ & -& \frac{2}{{(1+{d}_{j})}^{3}}(\frac{{\lambda}_{m}}{4}{(\mathrm{tr}(\epsilon )+{\varepsilon}_{zz})}^{2}{G}_{m}(\mathrm{tr}(\epsilon ),{d}_{j})+\frac{{\mu}_{m}}{2}(\sum_{i=1}^{2}{\tilde{\epsilon}}_{i}^{2}{G}_{m}({\tilde{\epsilon}}_{i},{d}_{j})))\\ & -& \frac{2\alpha (1-{\gamma}_{f})}{{(\alpha +{d}_{j})}^{3}}(\frac{{\lambda}_{f}}{2}\mathrm{tr}{(\kappa )}^{2}H({(-1)}^{(j+1)}\mathrm{tr}(\kappa ))+{\mu}_{f}\sum_{i=1}^{2}{\tilde{\kappa}}_{i}^{2}H({(-1)}^{(j+1)}{\tilde{\kappa}}_{i}))\end{array}\end{split}\]

où:

(3693)\[\frac{\partial {\varepsilon}_{zz}}{\partial {d}_{j}}=-\frac{{\lambda}_{m}(\mathrm{tr}(\epsilon )+{\varepsilon}_{zz})}{(2{\mu}_{m}+{\lambda}_{\mathrm{m.}}{\xi}_{m}(\mathrm{tr}(\epsilon ),{d}_{j}))}\frac{\partial {\xi}_{m}(\mathrm{tr}(\epsilon ),{d}_{j})}{\partial {d}_{j}}\]

et:

(3694)\[\frac{\partial {\xi}_{m}}{\partial {d}_{j}}(x,{d}_{j})=-\frac{1}{2}\left(\frac{1-{\gamma}_{\mathrm{mt}}}{{(1+{d}_{j})}^{2}}H(x)+\frac{{\alpha}_{c}(1-{\gamma}_{\mathrm{mc}})}{{({\alpha}_{c}+{d}_{j})}^{2}}H(-x) \right)=-\frac{{G}_{m}(x,{d}_{j})}{2{(1+{d}_{j})}^{2}}\]

On vérifie qu’on a bien \(\frac{\partial {\xi}_{m}}{\partial {d}_{j}}(x,{d}_{1,}{d}_{2})<0\) comme attendu.

Opérateur de raideur tangente

Comme l’objectif principal d’un modèle global est de proposer une approche simplifiée à la modélisation d’un matériau complexe, tel que le béton armé, il est essentiel que sa performance numérique soit optimale. Ainsi, pour rendre le modèle adapté aux calculs avec des schémas implicites en temps, soit en quasi-statique soit en dynamique transitoire, le calcul de la raideur tangente cohérente devient indispensable pour avoir une convergence quadratique et robuste du processus itératif global de Newton.

L’essentiel du calcul du modèle s’effectue dans le repère des vecteurs propres des tenseurs de déformations généralisées (et des contraintes généralisées, cf. (4899)), la raideur tangente étant donc elle aussi exprimée dans le même repère. La transformation nécessaire pour pouvoir ensuite l’utiliser dans l’assemblage de la matrice de raideur globale est précisée dans le chapitre suivant.

Afin de simplifier l’écriture on définit un vecteur de contraintes généralisées (effort membranaire, moment fléchissant), \(\Sigma\) , et un vecteur de déformations généralisées (extension, courbure), \(E\) , comme:

\(\Sigma ={({\tilde{N}}_{1}{\tilde{N}}_{2}{\tilde{M}}_{1}{\tilde{M}}_{2})}^{T}\)

\(E={({\tilde{\epsilon}}_{1}{\tilde{\epsilon}}_{2}{\tilde{\kappa}}_{1}{\tilde{\kappa}}_{2})}^{T}\)

L’opérateur de raideur tangente \(C\) est défini par la relation dans l’évolution réelle:

\(d\Sigma =C\cdot dE\)

Il peut être calculé comme la somme de deux contributions, celle qui correspond à une non-évolution de l’endommagement et celle qui est due à l’évolution de l’endommagement. Ces contributions peuvent être nommées: contribution élastique et contribution dissipative :

(3695)\[C=\underset{{C}_{e}}{\underset{\underbrace{}}{{\frac{d\Sigma }{dE}\mid }_{{\dot{d}}_{j}=0}}}+\underset{{C}_{d}}{\underset{\underbrace{}}{{\frac{d\Sigma }{dE}\mid }_{{f}_{{d}_{j}}=0}}}\]

De plus, on tient compte de la structure du tenseur \(C\) composé des contributions effort normal-extension, moment-courbure et de leurs couplages. Plus particulièrement les tenseurs \({C}_{e}\) et \({C}_{d}\) sont de la forme suivante :

(3696)\[\begin{split}{C}_{e}=\left(\begin{array}{cc}{C}_{e}^{\mathrm{mm}}& 0\\ 0& {C}_{e}^{\mathrm{ff}}\end{array}\right)\end{split}\]

On voit de (3696) que le couplage moments/extension et efforts de membrane/courbure n’est introduit qu’à travers la partie dissipative. Ce couplage a une justification physique, puisque toute fissuration perpendiculaire au feuillet moyen de la plaque affecte aussi bien le comportement en membrane qu’en flexion.

Les sous-matrices \({C}_{e}^{\mathrm{mm}}\) , \({C}_{e}^{\mathrm{ff}}\) , \({C}_{d}^{\mathrm{mm}}\) , \({C}_{d}^{\mathrm{mf}}\) , \({C}_{d}^{\mathrm{ff}}\) sont données dans le repère propre par les expressions qui suivent:

(3697)\[{({C}_{e}^{\mathrm{mm}})}_{ij}={\frac{\partial {\tilde{N}}_{i}}{\partial {\tilde{\epsilon}}_{j}}\mid }_{\dot{{d}_{k}}=0}=\frac{2{\lambda}_{m}{\xi}_{m}(\mathrm{tr}(\epsilon ),{d}_{k})}{2{\mu}_{m}+{\lambda}_{m}{\xi}_{m}(\mathrm{tr}(\epsilon ),{d}_{k})}+2{\mu}_{m}{\xi}_{m}({\tilde{\epsilon}}_{j},{d}_{k}){\delta}_{ij} {({C}_{e}^{\mathrm{ff}})}_{ij}={\frac{\partial {\tilde{M}}_{i}}{\partial {\tilde{\kappa}}_{j}}\mid }_{\dot{{d}_{k}}=0}={\lambda}_{f}{\xi}_{f}(\mathrm{tr}(\kappa ),{d}_{k})+2{\mu}_{f}{\xi}_{f}({\tilde{\kappa}}_{j},{d}_{k}){\delta}_{ij}\]

(3698)\[{({C}_{d}^{\mathrm{mm}})}_{ij}={\frac{\partial {\tilde{N}}_{i}}{\partial {\tilde{\epsilon}}_{j}}\mid }_{{f}_{{d}_{k}}=0}={\frac{\partial {\tilde{N}}_{i}}{\partial {d}_{k}}\cdot \frac{d({d}_{k})}{d{\tilde{\epsilon}}_{j}}\mid }_{{f}_{{d}_{k}}=0} {({C}_{d}^{\mathrm{ff}})}_{ij}={\frac{\partial {\tilde{M}}_{i}}{\partial {\tilde{\kappa}}_{j}}\mid }_{{f}_{{d}_{k}}=0}={\frac{\partial {\tilde{M}}_{i}}{\partial {d}_{k}}\cdot \frac{d({d}_{k})}{\partial {\tilde{\kappa}}_{j}}\mid }_{{f}_{{d}_{k}}=0} {({C}_{d}^{\mathrm{mf}})}_{ij}={\frac{\partial {\tilde{N}}_{i}}{\partial {\tilde{\kappa}}_{j}}\mid }_{{f}_{{d}_{k}}=0}={\frac{\partial {\tilde{N}}_{i}}{\partial {d}_{k}}\cdot \frac{d({d}_{k})}{\partial {\tilde{\kappa}}_{j}}\mid }_{{f}_{{d}_{k}}=0}\]

où

(3699)\[\frac{\partial {\tilde{N}}_{i}}{\partial {d}_{k}}={\lambda}_{m}(\mathrm{tr}(\epsilon )+{\varepsilon}_{zz})\frac{\partial {\xi}_{m}(\mathrm{tr}(\epsilon ),{d}_{j})}{\partial {d}_{k}}+{\lambda}_{m}{\xi}_{m}(\mathrm{tr}(\epsilon ),{d}_{j})\frac{\partial {\varepsilon}_{zz}}{\partial {d}_{k}}+2{\mu}_{m}{\tilde{\epsilon}}_{i}\frac{\partial {\xi}_{m}({\tilde{\epsilon}}_{i},{d}_{j})}{\partial {d}_{k}} \text{et} \quad \frac{\partial {\tilde{M}}_{i}}{\partial {d}_{k}}={\lambda}_{f}\mathrm{tr}(\kappa )\frac{\partial {\xi}_{f}(\mathrm{tr}(\kappa ),{d}_{1,}{d}_{2})}{\partial {d}_{k}}+2{\mu}_{f}{\tilde{\kappa}}_{i}\frac{\partial {\xi}_{f}({\tilde{\kappa}}_{i},{d}_{1,}{d}_{2})}{\partial {d}_{k}} \text{avec} \quad \frac{\partial {\xi}_{f}}{\partial {d}_{k}}(x,{d}_{1,}{d}_{2})=-\alpha \frac{(1-{\gamma}_{f})}{{(\alpha +{d}_{k})}^{2}}H({(-1)}^{(k+1)}x)\]

En plus des expressions (3699), on a recours aux équations (3694) pour déterminer \(\frac{\partial {\xi}_{m}}{\partial {d}_{k}}(x,{d}_{1,}{d}_{2})\) .

Les dérivées \({\frac{d({d}_{k})}{d{\tilde{\epsilon}}_{i}}\mid }_{{f}_{{d}_{k}}=0}\) et \({\frac{d({d}_{k})}{d{\tilde{\kappa}}_{i}}\mid }_{{f}_{{d}_{k}}=0}\) sont calculées en différentiant l’équation \({R}_{{d}_{k}}=0\) respectivement par rapport à \({d}_{k}\) , \({\tilde{\epsilon}}_{i}\) et \({\tilde{\kappa}}_{i}\) , (voir (3690)). Lorsque les deux mécanismes d’endommagement sont activé, on est amené à résoudre les systèmes qui suivent:

(3700)\[\begin{split}\left(\begin{array}{cc}{A}_{11}& {A}_{12}\\ {A}_{21}& {A}_{22}\end{array}\right)\left(\begin{array}{c}\frac{d({d}_{1})}{d{\tilde{\epsilon}}_{i}}\\ \frac{d({d}_{2})}{d{\tilde{\epsilon}}_{i}}\end{array}\right)=\left(\begin{array}{c}{B}_{1}^{m}\\ {B}_{2}^{m}\end{array}\right)\end{split}\]

avec \({A}_{\mathrm{jk}}=\frac{\partial {Y}_{j}}{\partial {d}_{k}}\) , \({B}_{k}^{m}=\frac{-\partial {Y}_{k}}{\partial \epsilon }=\frac{\partial \tilde{N}}{\partial {d}_{k}}\) , \({B}_{k}^{f}=\frac{-\partial {Y}_{k}}{\partial \kappa }=\frac{\partial \tilde{M}}{\partial {d}_{k}}\) ;

donc:

\(\frac{d({d}_{1})}{d{\tilde{\epsilon}}_{i}}=\frac{{A}_{22}{B}_{1}^{m}-{A}_{12}{B}_{2}^{m}}{{A}_{22}{A}_{11}-{A}_{12}{A}_{21}}\) et \(\frac{d({d}_{2})}{d{\tilde{\epsilon}}_{i}}=\frac{{A}_{11}{B}_{2}^{m}-{A}_{21}{B}_{1}^{m}}{{A}_{11}{A}_{22}-{A}_{21}{A}_{12}}\)

\(\frac{d({d}_{1})}{d{\tilde{\kappa}}_{i}}=\frac{{A}_{22}{B}_{1}^{f}-{A}_{12}{B}_{2}^{f}}{{A}_{22}{A}_{11}-{A}_{12}{A}_{21}}\) et \(\frac{d({d}_{2})}{d{\tilde{\kappa}}_{i}}=\frac{{A}_{11}{B}_{2}^{f}-{A}_{21}{B}_{1}^{f}}{{A}_{11}{A}_{22}-{A}_{21}{A}_{12}}\)

où:

(3701)\[\begin{split}\begin{array}{ccc}{A}_{\mathrm{jk}}& =& \frac{1}{{(1+{d}_{j})}^{2}}(\frac{{\lambda}_{m}}{2}(\mathrm{tr}(\epsilon )+{\varepsilon}_{zz})\frac{\partial {\varepsilon}_{zz}}{\partial {d}_{k}}{G}_{m}(\mathrm{tr}(\epsilon ),{d}_{j}))\\ & +& (\frac{1}{{(1+{d}_{j})}^{2}}(\frac{{\lambda}_{m}}{4}{(\mathrm{tr}(\epsilon )+{\varepsilon}_{zz})}^{2}\frac{\partial {G}_{m}(\mathrm{tr}(\epsilon ),{d}_{j})}{\partial {d}_{k}}+\frac{{\mu}_{m}}{2}\sum_{i=1}^{2}{\tilde{\epsilon}}_{i}^{2}\frac{\partial {G}_{m}({\tilde{\epsilon}}_{i},{d}_{j})}{\partial {d}_{k}})\\ & -& \frac{2}{{(1+{d}_{j})}^{3}}(\frac{{\lambda}_{m}}{4}{(\mathrm{tr}(\epsilon )+{\varepsilon}_{zz})}^{2}{G}_{m}(\mathrm{tr}(\epsilon ),{d}_{j})+\frac{{\mu}_{m}}{2}\sum_{i=1}^{2}{\tilde{\epsilon}}_{i}^{2}{G}_{m}({\tilde{\epsilon}}_{i},{d}_{j}))\\ & -& \frac{2\alpha (1-{\gamma}_{f})}{{(\alpha +{d}_{j})}^{3}}(\frac{{\lambda}_{f}}{2}\mathrm{tr}{(\kappa )}^{2}H({(-1)}^{(j+1)}\mathrm{tr}(\kappa ))+{\mu}_{f}\sum_{i=1}^{2}{\tilde{\kappa}}_{i}^{2}H({(-1)}^{(j+1)}{\tilde{\kappa}}_{i}))){\delta}_{\mathrm{jk}}\end{array}\end{split}\]

avec:

(3702)\[\frac{\partial {G}_{m}(x,{d}_{j})}{\partial {d}_{k}}=(\frac{2{\alpha}_{c}(1+{\gamma}_{\mathrm{mc}})(1+{d}_{j})({\alpha}_{c}-1)}{{({\alpha}_{c}+{d}_{j})}^{3}})H(-x)\]

et \(\frac{\partial {\varepsilon}_{zz}}{\partial {d}_{k}}\) donnés par (3690). La matrice \(({A}_{\mathrm{jk}})\) est bien inversible, cf. section numref:r7.01.32-integration-numerique .

Changement de repère

La démarche de changement de repère est identique à celle développée pour le modèle ENDO_ISOT_BETON (voir section 2.4.4.1 de [r7.01.04]) avec la seule différence qu’elle s’applique aux contraintes et déformations généralisées. Nous obtenons ainsi les composantes de:

\(C=(\begin{array}{cc}{C}^{\mathrm{mm}}& {C}^{\mathrm{mf}}\\ {({C}^{\mathrm{mf}})}^{T}& {C}^{\mathrm{ff}}\end{array})\)

comme étant:

\(\begin{array}{ccccc}{({C}^{\mathrm{mm}})}_{ijkl}& =& \frac{\partial {N}_{ij}}{\partial {\epsilon}_{kl}}& =& \sum_{m,n}{Q}_{\text{im}}^{m}{Q}_{\mathrm{jm}}^{m}{Q}_{\mathrm{kn}}^{m}{Q}_{\ln}^{m}\cdot \frac{\partial {\tilde{N}}_{m}}{\partial {\tilde{\epsilon}}_{n}}\\ & & & +& \frac{1}{2}\sum_{\begin{array}{}m,n\\ m\ne n\end{array}}(\frac{({Q}_{\mathrm{km}}^{m}{Q}_{\ln}^{m}+{Q}_{\mathrm{lm}}^{m}{Q}_{\mathrm{kn}}^{m})({Q}_{in}^{m}{Q}_{\mathrm{jm}}^{m}+{Q}_{\mathrm{jn}}^{m}{Q}_{\text{im}}^{m})}{{\tilde{\epsilon}}_{n}-{\tilde{\epsilon}}_{m}}){\tilde{N}}_{m}\end{array}\)

\(\begin{array}{ccccc}{({C}^{\mathrm{ff}})}_{ijkl}& =& \frac{\partial {M}_{ij}}{\partial {\kappa}_{kl}}& =& \sum_{m,n}{Q}_{\text{im}}^{f}{Q}_{\mathrm{jm}}^{f}{Q}_{\mathrm{kn}}^{f}{Q}_{\ln}^{f}\cdot \frac{\partial {\tilde{M}}_{m}}{\partial {\tilde{\kappa}}_{n}}\\ & & & +& \frac{1}{2}\sum_{\begin{array}{}m,n\\ m\ne n\end{array}}(\frac{({Q}_{\mathrm{km}}^{f}{Q}_{\ln}^{f}+{Q}_{\mathrm{lm}}^{f}{Q}_{\mathrm{kn}}^{f})({Q}_{in}^{f}{Q}_{\mathrm{jm}}^{f}+{Q}_{\mathrm{jn}}^{f}{Q}_{\text{im}}^{f})}{{\tilde{\kappa}}_{n}-{\tilde{\kappa}}_{m}}){\tilde{M}}_{m}\end{array}\)

\({({C}^{\mathrm{mf}})}_{ijkl}=\frac{\partial {N}_{ij}}{\partial {\kappa}_{kl}}=\sum_{m,n}{Q}_{\text{im}}^{m}{Q}_{\mathrm{jm}}^{m}{Q}_{\mathrm{kn}}^{f}{Q}_{\ln}^{f}\cdot \frac{\partial {\tilde{N}}_{m}}{\partial {\tilde{\kappa}}_{n}}\)

Les contraintes généralisées quant à elles s’écrivent comme, cf. (4899):

\(N={Q}_{m}\tilde{N}{Q}_{m}^{T}\)

\(M={Q}_{f}\tilde{M}{Q}_{f}^{T}\)

Calcul de la dissipation

Par définition, la densité de puissance de dissipation lors de l’endommagement vaut:

\(\dot{D}=\sigma :\dot{\varepsilon}-{\dot{\Phi}}_{\mathrm{ed}}^{S}=-\sum_{j=1,2}\frac{\partial {\Phi}_{\mathrm{ed}}^{S}}{\partial {d}_{j}}{\dot{d}}_{j}=\sum_{j=1,2}{Y}_{j}{\dot{d}}_{j}\)

Dans cette expression, on a utilisé la définition de \({Y}_{j}\) (3683). Dans la phase endommageante, les fonctions seuils satisfont toujours \({f}_{{d}_{j}}={Y}_{j}-{k}_{0}\equiv 0\) . Par conséquent, on peut calculer la dissipation cumulée comme:

(3703)\[D=\int\dot{D}\mathrm{dt}={k}_{0}({d}_{1}+{d}_{2})\]

On a démontré que la dissipation cumulée du processus d’endommagement est directement liée aux variables internes. Il suffit de faire la somme des deux contributions et de la multiplier par la constante de seuil \({k}_{0}\) .

Le calcul de la dissipation est réalisé par les options DISS_ELGA et DISS_ELNO de CALC_CHAMP. Ces champs ont une seule composante nommée ENDO.

Variables internes du modèle

Le modèle nécessite deux variables internes, \({d}_{1}\) et \({d}_{2}\) (correspondant aux variables \(\mathrm{V1}\) et \(\mathrm{V2}\) du Code_Aster ), qui représentent l’endommagement du côté de la face supérieure et du côté de la face inférieure, respectivement. La distinction entre les faces supérieure et inférieure s’effectue à travers l’orientation du repère local de chaque point de Gauss. Ainsi, la face inférieure s’endommage avec \(\dot{{d}_{2}}>0\) pour des courbures positives, et la face supérieure avec \(\dot{{d}_{1}}>0\) pour des courbures négatives, ce qui se déduit directement de la définition de \({\xi}_{f}\) , (cf. section r7.01.32-energie-libre).

En tout cas le choix de l’orientation du repère en un point de Gauss n’affecte pas le résultat final en ce qui concerne les déplacements et les rotations. Il peut y avoir un impact sur l’interprétation de l’endommagement et des contraintes généralisées si les orientations locales ne sont pas cohérentes dans une structure. Il est fortement conseillé d’assurer cette cohérence en renseignant le mot clé ANGL_REP avec les angles nautiques du repère local, (voir [u4.42.01]):

AFFE_CARA_ELEM(COQUE = _F(ANGL_REP = \((a,b)\) ))

En plus des variables \(\mathrm{V1}\) et \(\mathrm{V2}\) , on introduit aussi \(\mathrm{V3}\) et \(\mathrm{V4}\) , de valeur binaire (\(0\) ou \(1\) ), qui indiquent l’évolution instantanée de \(\mathrm{V1}\) et \(\mathrm{V2}\) . Plus précisément, \(\mathrm{V3}\) vaut \(1\) lorsque \(\mathrm{V1}\) évolue et \(0\) sinon. De même, \(\mathrm{V4}\) vaut 1 lorsque \(\mathrm{V2}\) évolue et \(0\) sinon.

On introduit aussi \(\mathit{V5}\) , \(\mathit{V6}\) et \(\mathit{V7}\) , dont le rôle est mesurer l’affaiblissement relatif de raideur de la dalle en béton armé de manière rationnelle, par exemple par visualisation en chaque point matériel:

\(\mathrm{V5}=1-\frac{1}{2}\left(\frac{1+{\gamma}_{\mathrm{mt}}{d}_{1}}{1+{d}_{1}}+\frac{1+{\gamma}_{\mathrm{mt}}{d}_{2}}{1+{d}_{2}}\right)\) \(\mathrm{V6}=1-\frac{1}{2}\left(\frac{{\alpha}_{c}+{\gamma}_{\mathrm{mc}}{d}_{1}}{{\alpha}_{c}+{d}_{1}}+\frac{{\alpha}_{c}+{\gamma}_{\mathrm{mc}}{d}_{2}}{{\alpha}_{c}+{d}_{2}}\right)\) et \(\mathrm{V7}=1-\mathrm{Max}\left(\frac{1+{\gamma}_{f}{d}_{1}}{1+{d}_{1}},\frac{1+{\gamma}_{f}{d}_{2}}{1+{d}_{2}}\right)\)

respectivement en traction, en compression et en flexion. Ces variables seront toujours comprises entre \(0\) et les \(1-{\gamma}_{.}\) respectifs, et seront toujours croissantes, étant nulles en absence d’endommagement. Ces variables sont plus «parlantes» que les variables \(\mathrm{V1}\) et \(\mathrm{V2}\) .

On introduit les variables internes \(\mathit{V8}\) , \(\mathit{V9}\) , \(\mathit{V10}\) et \(\mathit{V11}\) , afin de faciliter les post-traitements des armatures. Ces variables dépendent des déformations limites EPSI_ELS (limite à ELS) et EPSI_LIM (limite à l’ELU), qui sont fournies via DEFI_MATERIAU/GLRC_DM.

• \(V8\) : ACIXELS: rapport entre la déformation de l’acier dans la direction X (maximum entre la nappe inférieure et supérieure) et la déformation EPSI_ELS

• \(V9\) : ACIXELU: rapport entre la déformation de l’acier dans la direction X (maximum entre la nappe inférieure et supérieure) et la déformation EPSI_ELU. On rappelle que le modèle GLRC_DM ne modélise pas la plastification des aciers et n’est donc pas prédictif dans ce domaine.

• \(V10\) : ACIYELS: rapport entre la déformation de l’acier dans la direction Y (maximum entre la nappe inférieure et supérieure) et la déformation EPSI_ELS

• \(V11\) : ACIYELU: rapport entre la déformation de l’acier dans la direction Y(maximum entre la nappe inférieure et supérieure) et la déformation EPSI_ELU. On rappelle que le modèle GLRC_DMne modélise pas la plastification des aciers et n’est donc pas prédictif dans ce domaine.

Pour le béton, on introduit les variables internes \(\mathit{V12}\) et \(\mathit{V13}\) .

• \(V12\) : BETSUP: rapport entre la déformation principale la plus faible (en compression) du béton en face supérieure et la déformation limite du béton en compression EPSI_C.

• \(V13\) : BETINF: rapport entre la déformation principale la plus faible (en compression) du béton en face inférieure et la déformation limite du béton en compression EPSI_C.

Pour tracer les maximums des déformations, on introduit \(\mathit{V14}\) , \(\mathit{V15}\) et \(\mathit{V16}\) :

• \(V14\) : TRAMAX : déformation maximale temporelle en traction

• \(V15\) : COMMAX : déformation maximale temporelle en compression

• \(V16\) : FLEMAX: déformation maximale temporelleen flexion

Pour évaluer la validité des paramètres, on propose le calcul des erreurs suivantes:

• \(V17\) : ERRCOM: erreur exprimée en pourcentage entre les aires de la courbe en compression théorique (définie par EPSI_C et FCJ ) et la courbe approximée ( NYC, GAMMA_C ) entre 0 et COMMAX

• \(V18\) : ERRFLE: erreur exprimée en pourcentage entre les aires de la courbe en flexion théorique (calcul de section de béton armé)et la courbe approximée ( MFY, GAMMA_F) entre 0 et FLEMAX

Utilisation du modèle GLRC_DM en situations thermomécaniques

On admet que la distribution de température est affine dans l’épaisseur: cela est admissible pour une plaque en situation thermique stationnaire. Le coefficient de dilatation thermique est supposé être celui du béton, nettement majoritaire dans la section: il est donc directement utilisé dans les situations de membrane et de flexion de la plaque.

On admet aussi que les coefficients et paramètres matériaux ne dépendent pas de la température dans la gamme étudiée.

Cette distribution de température affine dans l’épaisseur se traduit donc par une déformation membranaire et une variation de courbure thermiques, d’où un simple décalage des déformations membranaires et de courbures en fonction des températures en paroi de la plaque.

Ainsi, une fois ce décalage opéré, il n’y a pas pas d’impact direct sur l’expression de la loi de comportement ni sur le calcul des évolutions irréversibles. Deux tests vérifient les résultats obtenus, cf. r7.01.32-verification .

Le modèle GLRC_DM peut donc être utilisé en situations thermomécaniques stationnaires sans modifications, dans ce cadre d’hypothèses.

Paramètres de la loi

Le modèle de dalle en béton armé endommageable GLRC_DM a donc besoin de paramètres caractéristiques de l’élasticité, complétés de 6 paramètres pour décrire le comportement d’endommagement: \({k}_{0}\) , pour définir la limite élastique, \(\alpha\) pour déterminer la participation de la flexion (voir r7.01.32-contraintes-generalisees ), \({\gamma}_{\mathrm{mt}}\) , \({\gamma}_{\mathrm{mc}}\) et \({\alpha}_{c}\) , \({\gamma}_{f}\) pour décrire la réponse non linéaire. Tous ces paramètres peuvent être identifiés à partir d’essais de traction et de flexion purs uniaxiaux monotones. Certains d’entre eux sont remplacés dans le jeu de données code_aster par des paramètres plus «parlants», voir ci-après au paragraphe r7.01.32-bilan-identification-parametres .

Il est possible de procéder soit à partir d’estimations analytiques simples (qui donnent les ordres de grandeur) soit à partir d’un recalage sur une courbe de réponse fournie par un autre modèle de comportement, éventuellement en intégrant des compromis.

On décrit dans les paragraphes ci-dessous la démarche et on en fait le bilan au paragraphe r7.01.32-bilan-identification-parametres .

Identification des paramètres de comportement élastique linéaire

Dans ce modèle on suppose que le milieu béton armé est homogénéisé et on laisse à l’utilisateur le soin de choisir (calculer ou mesurer) les paramètres: \({E}_{\mathrm{eq}}^{m}\) (module de Young effectif en membrane), \({E}_{\mathrm{eq}}^{f}\) (module de Young effectif en flexion), \({\nu}_{m}\) (coefficient de Poisson effectif en membrane) et \({\nu}_{f}\) (coefficient de Poisson effectif en flexion). On applique les relations suivantes pour déterminer les coefficients de Lamé \({\lambda}_{m}\) , \({\mu}_{m}\) et \({\lambda}_{f}\) , \({\mu}_{f}\) :

(3704)\[\begin{split}\begin{array}{ccc}{\lambda}_{m}=\frac{{\nu}_{m}{E}_{\mathrm{eq}}^{m}h}{(1+{\nu}_{m})(1-2{\nu}_{m})}& ,& {\mu}_{m}=\frac{{E}_{\mathrm{eq}}^{m}h}{2(1+{\nu}_{m})}\\ {\lambda}_{f}=\frac{{\nu}_{f}{E}_{\mathrm{eq}}^{f}{h}^{3}}{12(1-{\nu}_{f}^{2})}& ,& {\mu}_{f}=\frac{{E}_{\mathrm{eq}}^{f}{h}^{3}}{24(1+{\nu}_{f})}\end{array}\end{split}\]

Les relations ci-dessus ne sont pas interchangeables par F↔M pour les paramètres membrane et flexion, puisqu’en membrane les relations correspondent au cas général (élasticité 3D) et la condition de contraintes planes est traitée au sein de la formulation du modèle, tandis que pour la flexion on s’est placé d’emblée en élasticité 2D avec contraintes planes. Dans le domaine élastique on a ainsi:

\(\begin{array}{}{N}_{\alpha \beta }=\frac{{E}_{\mathrm{eq}}^{m}h}{1-{\nu}_{m}^{2}}({\nu}_{m}.\mathrm{tr}\epsilon .{\delta}_{\alpha \beta }+(1-{\nu}_{m}){\epsilon}_{\alpha \beta })\\ {M}_{\alpha \beta }=\frac{{E}_{\mathrm{eq}}^{f}{h}^{3}}{12(1-{\nu}_{f}^{2})}({\nu}_{f}.\mathrm{tr}\kappa .{\delta}_{\alpha \beta }+(1-{\nu}_{f}){\kappa}_{\alpha \beta })\end{array}\)

\(\alpha\) et \(\beta\) étant des indices allant de 1 à 2.

Par défaut \({E}_{\mathrm{eq}}^{m}={E}_{\mathrm{eq}}^{f}=E\) et \({\nu}_{m}={\nu}_{f}=\nu\) , où \(E\) et \(\nu\) sont les coefficients élastiques renseignés dans le fichier de commande sous le mot-clé ELAS. En revanche, comme le béton armé n’est pas un matériau homogène, la valeur effective de \({E}_{\mathrm{eq}}^{f}\) peut être différente de \({E}_{\mathrm{eq}}^{m}\) . Par conséquent, on laisse à l’utilisateur la possibilité d’introduire des valeurs \({E}_{\mathrm{eq}}^{f}\) et \({\nu}_{f}\) (EF et NUF sous le mot-clé facteur GLRC_DM) différentes de \(E\) et \(\nu\) , qui dans ce cas ne servent que pour décrire l’élasticité en membrane.

La condition des contraintes planes pour la membrane \({\sigma}_{zz}=0\) est satisfaite de la manière décrite dans le paragraphe r7.01.32-seuils-evolution-endommagement .

Remarque:

Dans (4900), on obtient une relation différente entre* \({\lambda}_{f}\) , \({\nu}_{f}\) et \({E}_{\mathrm{eq}}^{f}\) d’une part, et entre \({\lambda}_{m}\) , \({\nu}_{m}\) et \({E}_{\mathrm{eq}}^{m}\) d’autre part. Cette différence est directement liée à la prise encompte différente enmembrane et en flexionde la condition de contraintes planes. Plus particulièrement, on définit \({E}_{\mathrm{eq}}^{f}\) et \({\nu}_{f}\) à travers un essai de flexion pure, où \({\kappa}_{yy}=-{\nu}_{f}{\kappa}_{xx}\) , et \({M}_{ij}=0\) , sauf \({M}_{xx}\ne 0\) *. On se sert alors des équations suivantes pour trouver la relation entre* \({\lambda}_{f}\) , \({\mu}_{f}\) et \({E}_{\mathrm{eq}}^{f}\) , \({\nu}_{f}\) :

\({M}_{yy}=({\lambda}_{f}(1-{\nu}_{f})-2{\mu}_{f}{\nu}_{f}){\kappa}_{xx}=0\)

et

\({M}_{xx}=({\lambda}_{f}(1-{\nu}_{f})-2{\mu}_{f}){\kappa}_{xx}=\frac{{E}_{\mathrm{eq}}^{f}{h}^{3}}{12}{\kappa}_{xx}\)

d’où l’on obtient:

(3705)\[{E}_{\mathrm{eq}}^{f}=\frac{12}{{h}^{3}}({\lambda}_{f}(1-{\nu}_{f})+2{\mu}_{f}) \text{et} {\lambda}_{f}(1-{\nu}_{f})-2{\mu}_{f}{\nu}_{f}=0\]

En résolvant (4901) , on a les relations exprimées dans (4900).

L’identification des paramètres élastiques \({E}_{\mathrm{eq}}^{m}\) , \({\nu}_{m}\) , \({E}_{\mathrm{eq}}^{f}\) et \({\nu}_{f}\) du modèle à partir des caractéristiques du béton et des aciers repose sur deux cas de chargement: traction pure et flexion pure.

Considérons les caractéristiques suivantes pour le béton: module de Young \({E}_{b}\) , coefficient de Poisson \({\nu}_{b}\) , épaisseur de la dalle \(h\) , et pour les aciers: module de Young \({E}_{a}\) , coefficient de Poisson \({\nu}_{a}\) , section totale par mètre linéaire (pour les deux nappes, supposées symétriques dans l’épaisseur et identiques dans les deux directions) \({S}_{a}\) , position relative d’une nappe dans l’épaisseur \({\chi }_{a}\in \left]0,1\right[\) .

On obtient ainsi par l’essai uniaxial en traction élastique pure :

(3706)\[\begin{split}\begin{array}{}{N}_{xx}={E}_{\mathrm{eq}}^{m}h{\epsilon}_{xx}={E}_{a}{S}_{a}{\epsilon}_{xx}+\frac{{E}_{b}h}{1-{\nu}_{b}^{2}}({\epsilon}_{xx}+{\nu}_{b}{\epsilon}_{yy})\\ {N}_{yy}=0={E}_{a}{S}_{a}{\epsilon}_{yy}+\frac{{E}_{b}h}{1-{\nu}_{b}^{2}}({\epsilon}_{yy}+{\nu}_{b}{\epsilon}_{xx})\end{array}\end{split}\]

D’où (mots-clé E et NU):

(3707)\[{E}_{\mathrm{eq}}^{m}={E}_{a}\frac{{S}_{a}}{h}+{E}_{b}\cdot \frac{{E}_{b}h+{E}_{a}{S}_{a}}{{E}_{b}h+{E}_{a}{S}_{a}(1-{\nu}_{b}^{2})}\]

On observe que cette identification produit une erreur sur la raideur en cisaillement élastique plan de la dalle, cas pour lequel les aciers ne contribuent pas (il s’agit de grilles de tiges soudées), ce qui rend le comportement homogénéisé orthotrope et non pas isotrope. En effet, on obtient avec les valeurs (4903):

(3708)\[{G}_{\mathrm{eq}}^{m}=\frac{{E}_{\mathrm{eq}}^{m}}{2(1+{\nu}_{m})}=\frac{{E}_{b}}{2(1+{\nu}_{b})}\cdot \frac{{E}_{b}^{2}{h}^{2}+2{E}_{a}{E}_{b}h{S}_{a}+{E}_{a}^{2}{S}_{a}^{2}(1-{\nu}_{b}^{2})}{{E}_{b}^{2}{h}^{2}+{E}_{a}{E}_{b}h{S}_{a}(1-{\nu}_{b})}\ne \frac{{E}_{b}}{2(1+{\nu}_{b})}\]

Si l’on préfère assurer prioritairement l’identification sur le cas de cisaillement élastique plan de la dalle, et sur le cas de la réponse selon la direction de traction pure (donc en acceptant l’erreur sur l’effet de Poisson orthogonal), on obtient:

(3709)\[{E}_{\mathrm{eq}}^{m}={E}_{b}+\frac{{E}_{a}{S}_{a}(1-{\nu}_{b})}{h}\]

On prendra garde à ce que cette identification grossière (non admissible thermodynamiquement par rapport à l’essai de traction pure) ne donne pas des valeurs fantaisistes de \({\nu}_{m}\) .

Ensuite, on obtient par l’essai uniaxial en flexion élastique pure :

(3710)\[\begin{split}\begin{array}{}{M}_{xx}=\frac{{E}_{\mathrm{eq}}^{f}{h}^{3}}{12}{\kappa}_{xx}=\frac{1}{4}{E}_{a}{S}_{a}{h}^{2}{\chi }_{a}^{2}{\kappa}_{xx}+\frac{{E}_{b}{h}^{3}}{12(1-{\nu}_{b}^{2})}({\kappa}_{xx}+{\nu}_{b}{\kappa}_{yy})\\ {M}_{yy}=0=\frac{1}{4}{E}_{a}{S}_{a}{h}^{2}{\chi }_{a}^{2}{\kappa}_{yy}+\frac{{E}_{b}{h}^{3}}{12(1-{\nu}_{b}^{2})}({\kappa}_{yy}+{\nu}_{b}{\kappa}_{xx})\end{array}\end{split}\]

D’où (mots-clé EF et NUF):

(3711)\[{E}_{\mathrm{eq}}^{f}=\frac{3}{h}{E}_{a}{S}_{a}{\chi }_{a}^{2}+{E}_{b}\cdot \frac{{E}_{b}h+3{E}_{a}{S}_{a}{\chi }_{a}^{2}}{{E}_{b}h+3{E}_{a}{S}_{a}{\chi }_{a}^{2}(1-{\nu}_{b}^{2})}\]

On observe aussi que cette identification produit une erreur sur la raideur en flexion élastique anticlastique \({M}_{xy}\) de la dalle (coefficient \({G}_{\mathrm{eq}}^{f}=\frac{{E}_{\mathrm{eq}}^{f}{h}^{3}}{24(1+{\nu}_{f})}\) au lieu de \({G}_{b}^{f}=\frac{{E}_{b}{h}^{3}}{24(1+{\nu}_{b})}\) ), cas pour lequel les aciers ne contribuent pas.

Identification des paramètres de comportement élastique endommageable

La façon dont on obtient les paramètres de l’élasticité linéaire étant présentée au paragraphe r7.01.32-parametres-elastique-lineaire , on propose de calculer les paramètres d’endommagement du modèle à partir de trois tests: un essai de traction pure, un essai de compression pure et un essai de flexion pure monotones uniaxiaux.

De cette manière on obtient les valeurs du seuil \({k}_{0}\) , des trois paramètres relatifs aux effets de membrane (\({\gamma}_{\mathrm{mt}}\) , \({\alpha}_{c}\) et \({\gamma}_{\mathit{mc}}\) ) indépendamment des deux paramètres relatifs aux effets de flexion (\(\alpha\) , \({\gamma}_{f}\) ).

Paramètres de traction (mots-clé NYT,GAMMA_T)

Seuil d’apparition de l’endommagement en traction NYT

En particulier, pour la traction élastique uniaxiale pure à l’apparition de l’endommagement on peut écrire la valeur du seuil, cf. (3683):

\({f}_{{d}_{j}}={Y}_{j}^{m}-{k}_{0}={\epsilon}_{D}^{2}(\frac{{\lambda}_{m}}{4}{(1-2{\nu}_{m})}^{2}(1-{\gamma}_{\mathrm{mt}})+\frac{{\mu}_{m}}{2}(1-{\gamma}_{\mathrm{mt}}+{\nu}_{m}^{2}\frac{1-{\gamma}_{\mathrm{mc}}}{{\alpha}_{c}}))-{k}_{0}=0\)

\({\epsilon}_{D}\) étant la déformation élastique à l’apparition de l’endommagement, ayant alors \({\epsilon}_{yy}=-{\nu}_{m}{\epsilon}_{D}={\varepsilon}_{zz}\) et \({\xi}_{m}(x,0,0)=1\) , d’où:

(3712)\[{k}_{0}=\frac{{\lambda}_{m}{(1-2{\nu}_{m})}^{2}(1-{\gamma}_{\mathrm{mt}})+2{\mu}_{m}(1-{\gamma}_{\mathrm{mt}}+{\nu}_{m}^{2}(1-{\gamma}_{\mathrm{mc}})/{\alpha}_{c})}{4{({\lambda}_{m}(1-2{\nu}_{m})+2{\mu}_{m})}^{2}}{N}_{D}^{2} {k}_{0}=\frac{{N}_{D}^{2}}{4{E}_{\mathrm{éq}}^{m}h(1+{\nu}_{m})}\cdot ((1-{\nu}_{m})(1+2{\nu}_{m})(1-{\gamma}_{\mathrm{mt}})+{\nu}_{m}^{2}(1-{\gamma}_{\mathrm{mc}})/{\alpha}_{c})\]

ayant

\({N}_{D}=({\lambda}_{m}(1-2{\nu}_{m})+2{\mu}_{m}){\epsilon}_{D}={E}_{\mathrm{eq}}^{m}h{\epsilon}_{D}\)

Remarque

On rappelle que \({\gamma}_{\mathrm{mt}}\le 1\) (cf. (4289)) et qu* \({\gamma}_{\mathrm{mc}}\le 1\) pour que l’endommagement se traduise bien par un affaiblissement de la raideur. On observe aussi sur (4288) que l’on ne peut pas avoir à la fois \({\gamma}_{\mathrm{mt}}=1\) et \({\gamma}_{\mathrm{mc}}=1\) , car alors \({k}_{0}=0\) (le modèle n’a pas de domaine élastique), ou alors il faudrait donner \({N}_{D}=\infty\) (mot-clé NYT).

En poursuivant l’analyse faite au r7.01.32-parametres-elastique-lineaire , juste à l’apparition de l’endommagement, la contrainte longitudinale dans le béton vaut:

\({E}_{b}{\epsilon}_{D}\frac{1-{\nu}_{b}{\nu}_{m}}{1-{\nu}_{b}^{2}}\)

de telle sorte que l’on puisse exprimer le seuil \({N}_{D}\) (mot-clé NYT) avec la limite de fissuration du béton \({\sigma}_{b}^{t}\) en traction, en supposant valide le critère local \({\sigma}_{xx}\le {\sigma}_{b}^{t}\) :

(3713)\[{N}_{D}={\sigma}_{b}^{t}\frac{{E}_{\mathrm{eq}}^{m}h}{{E}_{b}}\cdot \frac{1-{\nu}_{b}^{2}}{1-{\nu}_{b}{\nu}_{m}}\]

Paramètre de dégradation de raideur en traction GAMMA_T

Le paramètres \({\gamma}_{\mathit{mt}}\) est égal au rapport entre la pente correspondant à la rigidité de la phase endommagée, fournie par le mot-clé PENTE/TRACTIONdans DEFI_GLRC et la pente correspondant à la rigidité élastique.

Pour rappel (voir u4.42.06), trois méthodes dénommées RIGI_ACIER, PLAS_ACIER et UTIL sont disponibles. Ces trois calculs de pentes permettent de mettre en place trois méthodes de recalage différentes en fonction des propriétés matériaux renseignées pour la traction. Dans le cas où la limite élastique des aciers n’est pas connue, les méthodes de recalage RIGI_ACIER, i.e. pente post-élastique égale à la pente de reprise de raideur des aciers, et UTIL, i.e. pente post-élastique coupe la pente de reprise de raideur des aciers à une déformation maximale dont la valeur est imposée par l’utilisateur, sont accessibles (cette méthode n’est pas adaptée pour des déformations maximales plus faibles que le creux de la courbe de référence, voir figure Fig. 253). Dans le cas où la limite d’élasticité des aciers est connue, il est possible d’utiliser la méthode de recalage à la limite de plasticité des aciers (PLAS_ACIER). Les différentes méthodes de recalage sont illustrées par les figures qui suivent.

Fig. 252 Courbe de traction (GLRC_DM vs Référence) Recalage PENTE = RIGI_ACIER

Fig. 253 Courbe de traction (GLRC_DM vs Référence) Recalage PENTE = UTIL

Fig. 254 Courbe de traction (GLRC_DM vs Référence) Recalage PENTE = PLAS_ACIER

Cas distorsion uniaxiale pure

Vérifions l’effet de ces paramètres sur l’apparition de l’endommagement suite à un chargement de distorsion pure élastique \({\epsilon}_{xy}=\tilde{{\epsilon}_{1}}=-\tilde{{\epsilon}_{2}}\) , avec \({\epsilon}_{xx}={\epsilon}_{yy}=0\) . Ainsi: \({N}_{xy}=\tilde{{N}_{1}}=-\tilde{{N}_{2}}=2{\mu}_{m}{\epsilon}_{xy}\) . Le seuil d’apparition de l’endommagement (3683) avec le modèle GLRC_DM est atteint pour l’effort de cisaillement:

\({N}_{xy}^{D}=2\frac{\sqrt{2{\mu}_{m}{k}_{0}}}{\sqrt{1+{\alpha}_{c}-{\gamma}_{\mathrm{mc}}-{\gamma}_{\mathrm{mt}}{\alpha}_{c}}}=\frac{{N}_{D}}{1+{\nu}_{m}}\cdot \sqrt{\frac{(1-{\nu}_{m})(1+2{\nu}_{m})(1-{\gamma}_{\mathrm{mt}})+{\nu}_{m}^{2}(1-{\gamma}_{\mathrm{mc}})/{\alpha}_{c}}{1+{\alpha}_{c}-{\gamma}_{\mathrm{mc}}-{\gamma}_{\mathrm{mt}}{\alpha}_{c}}}\)

Il peut être utile de confronter cette prédiction avec celle que l’on obtient avec le modèle de béton ENDO_ISOT_BETON [r7.01.04], avec la limite de fissuration du béton \({\sigma}_{b}^{t}\) en traction:

\({{N}_{xy}^{D}}_{\mathrm{EIB}}=2{\sigma}_{b}^{t}h\sqrt{\frac{(1-{\nu}_{b})(1+2{\nu}_{b})}{{(1+{\nu}_{b})}^{2}}}=2{N}_{D}\frac{{E}_{b}(1-{\nu}_{b}{\nu}_{m})}{{E}_{\mathrm{éq}}^{m}(1-{\nu}_{b}^{2})}\cdot \sqrt{\frac{(1-{\nu}_{b})(1+2{\nu}_{b})}{{(1+{\nu}_{b})}^{2}}}\)

Remarque

Dans toute situation combinée (compression+cisaillement, etc.) en membranaire pur, l’expression du seuil de premier endommagement \({Y}_{j}^{m}={k}_{0}\) du modèle GLRC_DM, cf. (3683), en fonction des efforts membranaires \({N}_{xx},{N}_{xy}\mathrm{...}\) , est constituée des mêmes monômes que le critère «habituel» en contraintes planes du matériau béton considéré ne pas résister au-delà de \({\sigma}_{b}^{t}\) . Ceci provient du choix de la formulation du modèle GLRC_DMen membrane, en filiation directe du modèle ENDO_ISOT_BETON.

Paramètres de compression (mots-clé NYC,GAMMA_C, ALPHA_C)

En compression, trois paramètres sont à déterminer: \({N}_{C}\) , \({\alpha}_{c}\) et \({\gamma}_{\mathrm{mc}}\) . On utilise trois équations pour y parvenir.

La première nous vient de la formulation du modèle. Si on envisage un essai de compression uniaxiale pure , la valeur du seuil de premier endommagement, cf. (3683), s’écrit en désignant par \({N}_{C}\) l’effort normal correspondant:

\({k}_{0}=\frac{{N}_{C}^{2}}{4{E}_{\mathrm{éq}}^{m}h(1+{\nu}_{m})}\cdot ((1-{\nu}_{m})(1+2{\nu}_{m})\frac{1-{\gamma}_{\mathrm{mc}}}{{\alpha}_{c}}+{\nu}_{m}^{2}(1-{\gamma}_{\mathrm{mt}}))\)

On doit donc nécessairement avoir la relation:

\(\frac{{N}_{C}^{2}}{{N}_{D}^{2}}=\frac{{\alpha}_{c}(1-{\nu}_{m})(1+2{\nu}_{m})(1-{\gamma}_{\mathrm{mt}})+{\nu}_{m}^{2}(1-{\gamma}_{\mathrm{mc}})}{(1-{\nu}_{m})(1+2{\nu}_{m})(1-{\gamma}_{\mathrm{mc}})+{\alpha}_{c}{\nu}_{m}^{2}(1-{\gamma}_{\mathrm{mt}})}\)

Ainsi, on obtient l’expression de \({\alpha}_{c}\) en fonction de \({\gamma}_{\mathrm{mc}}\) , \({\gamma}_{\mathrm{mt}}\) , et \({N}_{D}\) , \({N}_{C}\) :

(3714)\[{\gamma}_{\mathrm{mc}}=1-{\alpha}_{c}(1-{\gamma}_{\mathrm{mt}})\cdot \frac{{N}_{D}^{2}(1-{\nu}_{m})(1+2{\nu}_{m})-{N}_{C}^{2}{\nu}_{m}^{2}}{{N}_{C}^{2}(1-{\nu}_{m})(1+2{\nu}_{m})-{N}_{D}^{2}{\nu}_{m}^{2}} \iff {\alpha}_{c}=\frac{(1-{\gamma}_{\mathrm{mc}})}{(1-{\gamma}_{\mathrm{mt}})}\cdot \frac{{N}_{C}^{2}(1-{\nu}_{m})(1+2{\nu}_{m})-{N}_{D}^{2}{\nu}_{m}^{2}}{{N}_{D}^{2}(1-{\nu}_{m})(1+2{\nu}_{m})-{N}_{C}^{2}{\nu}_{m}^{2}}\]

Remarque:

Il est nécessaire que \({\gamma}_{\mathrm{mc}}\le 1\) , tout comme \({\gamma}_{\mathrm{mt}}\le 1\) , cf. (3683) . On rappelle aussi, cf. (4288) , que l’on ne peut pas avoir à la fois \({\gamma}_{\mathrm{mt}}=1\) et \({\gamma}_{\mathrm{mc}}=1\) . De l” (4290) , on obtient la condition nécessaire, quel que soit \({\alpha}_{c}\) :

\(\mid {N}_{C}\mid \le {N}_{D}\frac{\sqrt{(1-{\nu}_{m})(1+2{\nu}_{m})}}{{\nu}_{m}}\)

L’égalité dans la relation ci-dessus conduit à \({\gamma}_{\mathrm{mc}}=1\) .

Pour du béton armé, ayant \({\nu}_{m}\approx 0,2\) , cette condition s’écrit: \(\mid {N}_{C}\mid <5,2{N}_{D}\) .

La deuxième équation provient d’une condition supplémentaire qu’on impose: on vise le même endommagement maximal en traction et en compression. 𝑑𝑚𝑎𝑥 correspond donc à la perte de raideur la plus importante, à la fois en traction(noté 𝜉𝑡_𝑚𝑖𝑛) et en compression(noté 𝜉𝑐_𝑚𝑖𝑛).

Enfin, on introduit une variable \({N}_{\mathit{CU}}\) qui correspond à la résistance ultime de la section en compression, à partir de la hauteur h de la coque, et de la valeur de la résistance moyenne en compression \({F}_{\mathit{CJ}}\) et \({ϵ}_{\mathit{C1}}\) la déformation correspondante.

\({N}_{\mathit{CU}}=h{F}_{\mathit{CJ}}+{E}_{a}{S}_{a}\cdot {ϵ}_{\mathit{C1}}\)

(3715)\[\]

Le seuil de compression \({N}_{C}\) est alors lié à la limite ultime de compression \({N}_{\mathit{CU}}\) via l’équation:

(3716)\[{N}_{C}=\frac{{N}_{\mathit{CU}}-{\gamma}_{\mathit{mc}}{\mathit{hE}}_{m}{ϵ}_{\mathit{C1}}}{1-{\gamma}_{\mathit{mc}}}\]

On vise le même endommagement maximal en traction et en compression. \({d}_{max}\) doit donc permettre d’atteindre la perte de raideur la plus importante, à la fois en traction (noté \({\xi}_{{t}_{min}}\) ) et en compression (noté \({\xi}_{{c}_{min}}\) ) . En compression, l’endommagement maximal permet d’atteindre le pic de compression \({N}_{\mathit{CU}}\) . En traction, il permet d’atteindre le point de plastification des aciers.

(3717)\[\frac{{\alpha}_{c}+{\gamma}_{\mathit{mc}}{d}_{max}}{{\alpha}_{c}+{d}_{max}}={\xi}_{{c}_{min}}=\frac{{N}_{\mathit{CU}}}{{\mathit{hE}}_{\mathit{eq}}^{m}{ϵ}_{\mathit{C1}}}\]

(3718)\[\frac{1+{\gamma}_{\mathit{mt}}{d}_{max}}{1+{d}_{max}}={\xi}_{{t}_{min}}=\frac{{S}_{a}{\sigma}_{\mathit{ya}}}{{\mathit{hE}}_{m}{ϵ}_{\mathit{tmax}}}=\frac{{S}_{a}{\sigma}_{\mathit{ya}}}{{N}_{D}+({S}_{a}{\sigma}_{\mathit{ya}}-{N}_{D})/{\gamma}_{\mathit{mt}}}\]

La résolution de l’équation (4294) donne la valeur de l’endommagement maximal \({d}_{max}\) .

\({d}_{max}=\frac{{S}_{a}{\sigma}_{\mathit{ya}}-{N}_{D}}{{N}_{D}\ast {\gamma}_{\mathit{mt}}}\)

Dans le cas où la valeur obtenue pour \({d}_{max}\) est positive, on obtient:

\({\alpha}_{c}=\frac{{d}_{max}{N}_{C}(1-{\gamma}_{\mathit{mc}})}{h{E}_{m}{ϵ}_{\mathit{C1}}-{N}_{\mathit{CU}}}\) if faut avoir nécessairement \({N}_{\mathit{CU}}<h{E}_{m}{ϵ}_{\mathit{C1}}\)

et

\(\frac{{d}_{max}{N}_{C}}{h{E}_{m}{ϵ}_{\mathit{C1}}-{N}_{\mathit{CU}}}=\frac{({N}_{D}^{2}{\nu}_{m}^{2}-{N}_{C}^{2}(1-{\nu}_{m})(1+2{\nu}_{m}))}{(1-{\gamma}_{\mathit{mt}})({N}_{C}^{2}{\nu}_{m}^{2}-{N}_{D}^{2}(1-{\nu}_{m})(1+2{\nu}_{m}))}\)

Cette dernière équation est une équation cubique \({N}_{C}\)

Dans le cas où \({d}_{max}<0\) , on modifie le calcul de \({\gamma}_{\mathit{mt}}\) , en imposant \({\gamma}_{\mathit{mt}}=0\)

et l’équation (4294) est modifiée par

\(\frac{1}{1+{d}_{max}}=\frac{{S}_{a}{\sigma}_{\mathit{ya}}}{{\mathit{hE}}_{m}{ϵ}_{\mathit{tmax}}}\)

avec \({ϵ}_{\mathit{tmax}}=\frac{{\sigma}_{\mathit{ya}}}{{E}_{a}}\)

alors on a \({d}_{max}=\frac{{\mathit{hE}}_{m}}{{S}_{a}{\sigma}_{\mathit{ya}}}-1\) .

Paramètres de flexion (mots-clé MYF,GAMMA_F)

Options RIGI_ACIER et PLAS_ACIER

En flexion élastique uniaxiale pure un seul mécanisme d’endommagement est activé, selon son sens, positif ou négatif. Ici on choisit la flexion positive, pour laquelle on a toujours \({f}_{{d}_{1}}>{f}_{{d}_{2}}\) . La valeur maximale de courbure élastique \({\kappa}_{xx}\) à l’apparition de l’endommagement est notée \({\kappa}_{D}\) (\({\kappa}_{yy}=-{\nu}_{f}{\kappa}_{xx}\) ), telle que seul le seuil \({f}_{{d}_{1}}=0\) peut être atteint, tandis que \({f}_{{d}_{2}}<0\) pour chaque point de cette trajectoire du chargement, cf. (3683):

\({f}_{{d}_{1}}={Y}_{1}^{f}-{k}_{0}={\kappa}_{D}^{2}\frac{(1-{\gamma}_{f})}{\alpha}(\frac{{\lambda}_{f}}{2}{(1-{\nu}_{f})}^{2}+{\mu}_{f})-{k}_{0}\)

d’où:

(3719)\[\alpha =(1-{\gamma}_{f})\frac{{\lambda}_{f}{(1-{\nu}_{f})}^{2}+2{\mu}_{f}}{2{({\lambda}_{f}(1-{\nu}_{f})+2{\mu}_{f})}^{2}}\frac{{M}_{D}^{2}}{{k}_{0}}\]

ayant

\({M}_{D}=({\lambda}_{f}(1-{\nu}_{f})+2{\mu}_{f}).{\kappa}_{D}=\frac{{E}_{\mathrm{eq}}^{f}{h}^{3}}{12}{\kappa}_{D}\)

Comme la plaque en béton armé est supposée symétrique par rapport au feuillet moyen, on n’a besoin de faire l’identification que pour la flexion positive (la flexion négative donnant la même valeur).

En poursuivant l’analyse faite au paragraphe r7.01.32-parametres-elastique-lineaire , juste à l’apparition de l’endommagement, la contrainte longitudinale dans le béton vaut en paroi de la plaque (on sait qu’alors l’endommagement progresse immédiatement dans une bonne part de l’épaisseur de la section):

\({E}_{b}{\kappa}_{D}h\frac{1-{\nu}_{b}{\nu}_{f}}{2(1-{\nu}_{b}^{2})}\)

de telle sorte que l’on puisse exprimer le seuil \({M}_{D}\) avec la limite de fissuration du béton \({\sigma}_{t}^{b}\) :

(3720)\[{M}_{D}={\sigma}_{t}^{b}\frac{{E}_{\mathrm{eq}}^{f}{h}^{2}}{6{E}_{b}}\cdot \frac{1-{\nu}_{b}^{2}}{1-{\nu}_{b}{\nu}_{f}}\]

Cette valeur est utilisée pour les paramétrages RIGI_ACIER et PLAS_ACIER de de l’opérande PENTE/FLEXION de DEFI_GLRC.

Le paramètres \({\gamma}_{f}\) est égal au rapport entre la pente correspondant à la rigidité de la phase endommagée, fournie par le mot-clé PENTE/FLEXIONdans DEFI_GLRC et la pente correspondant à la rigidité élastique.Les options RIGI_ACIER et PLAS_ACIER correspondent au comportement présenté pour la traction au paragraphe r7.01.32-parametre-raideur-traction .

Afin de représenter plus fidèlement le comportement du béton armé en flexion pour une gamme de courbure plus restreinte,deux méthodes sont introduites pour l’identification des paramètres en flexion.

Ces deux méthodes sont basées sur l’identification au préalable de la courbe théorique courbure-flexion de la section. Cette courbe théorique est déterminée à partir de la théorique de la section plane.

La réponse de la section est calculée en intégrant les réponses sur la hauteur de la section.

Les hypothèses de réponse pour chaque matériau sont les suivantes:

• Réponse non-linéaire du béton définie par la limite en traction SYT, la limite en compression SYC et le paramètre c qui détermine la valeur de la pente du comportement en traction du béton après fissuration. On retient c=5, c’est à dire que la pente après fissuration vaut -0,2 \({p}_{\mathit{elas}}\) , où \({p}_{\mathit{elas}}\) est la pente élastique. Il s’agit d’un compromis évalué à partir d’essais numériques de traction avec la loi ENDO_ISOT_BETON pour des classes de béton couramment utilisées par l’ingénierie. On ajoute l’hypothèse que la limite en compression est supérieure à la limite en traction (en valeur absolue) \(\mid \mathit{SYC}\mid >\mid \mathit{SYT}\mid\) .

Fig. 255 Comportement du béton pour le calcul de la courbe courbure-moment

• Réponse élasto-plastique de l’acier défini par le paramètre SY

Fig. 256 Comportement du béton pour le calcul de la courbe courbure-moment

Pour déterminer la réponse courbure-flexion, on part de l’hypothèse d’une résultante d’effort axial nul. Pour la donnée d’une courbure, on détermine la valeur de la déformation à l’axe neutre afin d’obtenir une résultante d’effort axial nul. Cette détermination est réalisée en faisant des hypothèse sur les états possibles (fissuration du béton, limite en compression atteinte, plastification des aciers …etc). Elle ne fait pas appel à une intégration numérique avec un critère de précision, ce qui rend son utilisation plus sûre et moins coûteuse.

Le moment est ensuite déterminé à partir des valeurs de la déformation à l’axe neutre et la courbure.

Remarque:

Dans le cas de la modélisation de la section de la coque proposée, le module de flexion obtenu vaut:

\({E}_{\mathit{eq}}^{f}=\frac{3}{h}{E}_{a}{S}_{a}{\chi }_{a}^{2}+{E}_{b}\)

Il diffère donc du module qui prend en compte le coefficient de Poisson (4287).

Pour assurer la cohérence des formules et d’obtenir les même résultats en flexion élastique , on utilise un module de Young du béton modifié dans le calcul de la réponse de la section via l’intégration dans la hauteur.

\(\stackrel{̃}{{E}_{b}}={E}_{b}\cdot \frac{{E}_{b}h+3{E}_{a}{S}_{a}{\chi }_{a}^{2}}{{E}_{b}h+3{E}_{a}{S}_{a}{\chi }_{a}^{2}(1-{\nu}_{b}^{2})}\)

Deux options sont disponibles pour la détermination de l’approximation bilinéaire de courbe courbe moment-courbure .

Option RIGI_INIT

Pour l’option d’identification RIGI_INIT,on définit le moment de flexion seuil correspondant à l’apparition de l’endommagement de la façon suivante:

• A partir du point d’initiation de la première fissure \({M}_{D}\) , on recherche le plus petit moment sur la courbe élastique, pour lequel l’écart à la courbe théorique dépasse 5%. Ce moment est retenu pour la valeur \(\mathit{MYF}\) .

• La pente d’endommagement en flexion pure uniaxiale est alors définie comme la pente de la tangente à la courbe théorique moment-courbe et passant par le point seuil défini.

• Le paramètre \({\gamma}_{f}\) est alors défini comme le rapport des pentes endommageante et élastique.

Remarque:

Pour des sections avec peu de ferraillage, on peut obtenir une pente endommagente négative. \({\gamma}_{f}<0\) . Dans ce cas un message d’alarme est affiché, on impose une pente \({\gamma}_{f}={10}^{-3}\) et le moment seuil \(\mathit{MYF}\) est réduit pour satisfaire la condition de tangente à la courbe théorique.

Fig. 257 Approximation bilinéaire en flexion avec l’option RIGI_INIT

Pour l’implémentation de cette option, les étapes suivantes sont réalisées:

• Pour chaque courbure, dans un intervalle de courbures \(\left[{\kappa}_{D},{\kappa}_{L}\right]\) , on détermine le moment théorique \(M\) et la tangente \(\mathit{dM}/d\kappa\) à la courbe (\(d\kappa =({\kappa}_{L}-{\kappa}_{D})/1000\) ). Cette tangente est obtenue par différence finie.

• Parmi les valeurs calculées, on sélectionne la première courbure pour lequel l’écart à la courbe est supérieur à 5%.

• Pour obtenir la pente d’endommagement, on sélectionne la courbure limite \({\kappa}_{R}\) pour laquelle on minimise la différence entre \(\mathit{dM}/d\kappa\) et la pente d’endommagement.

Option UTIL

Pour l’option d’identificationUTIL, l’utilisateur définit une valeur de courbure cible via le mot-cléKAPPA_FLEX. On détermine alors la valeur de la courbure et du moment sueil \({\kappa}_{s}\) et \(\mathit{MYF}\) de façon à minimiser l’aire entre la courbe théorique et la courbe bilinéaire.

Remarque

Pour des sections avec peu de ferraillage, on peut obtenir une pente endommagement négative. \({\gamma}_{f}<0\) . Dans ce cas un message d’alarme est affiché, on impose une pente \({\gamma}_{f}={10}^{-3}\) et le moment seuil \(\mathit{MYF}\) est défini à la limite de fissuration du béton, i.e. \({M}_{D}\) .

Fig. 258 Approximation bilinéaire en flexion avec l’option UTIL

Pour l’implémentation de cette option, les étapes suivantes sont réalisées:

• Pour chaque courbure, dans un intervalle de courbures \(\left[{\kappa}_{D},{\kappa}_{\mathit{flex}}\right]\) (1000 points), on détermine le moment théorique \(M\) et on calcul l’aire sous la courbe théorique par la méthode des rectangles, ainsi que l’aire sous la courbe bilinéaire.

• On évalue l’aire sous la courbe des approximations bilinéaires qui garantissent la pente élastique calculée préalablement et qui coupent la courbe théorique pour \(\kappa ={\kappa}_{\mathit{flex}}\) : le point seuil est déplacé le long de la pente élastique et la pente endommagée ajustée pour passer par la courbe théorique.

• On retient l’approximation bilinéaire qui minimise l’écart entre l’aire sous sa courbe et celle de la courbe théorique.

Cas traction-flexion uniaxiales

Vérifions l’effet de ces paramètres sur l’apparition de l’endommagement suite à un chargement mêlant traction monotone uniaxiale et flexion monotone uniaxiale concomitants. Les seuils s’écrivent, pour chaque variable d’endommagement:

\(\begin{array}{}{k}_{0}\ge {Y}_{1}=\\ {\epsilon}_{xx}^{2}(\frac{{\lambda}_{m}}{4}{(1-2{\nu}_{m})}^{2}(1-{\gamma}_{\mathrm{mt}})+\frac{{\mu}_{m}}{2}(1-{\gamma}_{\mathrm{mt}}+{\nu}_{m}^{2}\frac{1-{\gamma}_{\mathrm{mc}}}{{\alpha}_{c}}))+{\kappa}_{xx}^{2}{\nu}_{f}^{2}{\mu}_{f}\frac{(1-{\gamma}_{f})}{\alpha}\end{array}\)

\(\begin{array}{}{k}_{0}\ge {Y}_{2}=\\ {\epsilon}_{xx}^{2}(\frac{{\lambda}_{m}}{4}{(1-2{\nu}_{m})}^{2}(1-{\gamma}_{\mathrm{mt}})+\frac{{\mu}_{m}}{2}(1-{\gamma}_{\mathrm{mt}}+{\nu}_{m}^{2}\frac{1-{\gamma}_{\mathrm{mc}}}{{\alpha}_{c}}))+{\kappa}_{xx}^{2}\frac{(1-{\gamma}_{f})}{\alpha}(\frac{{\lambda}_{f}}{2}{(1-{\nu}_{f})}^{2}+{\mu}_{f})\end{array}\)

On vérifie aisément que \({Y}_{1}<{Y}_{2}\) : le premier endommagement apparaît sur la variable \({d}_{2}\) comme attendu. Exploitons les résultats (4288) et (4290), le seuil \({Y}_{2}={k}_{0}\) s’écrit donc:

(3721)\[\frac{{N}_{xx}^{2}}{{N}_{D}^{2}}+\frac{{M}_{xx}^{2}}{{M}_{D}^{2}}=1\]

Il définit ainsi le domaine élastique prédit par le modèle GLRC_DM dans le quadrant \(({N}_{xx},{M}_{xx})\) positifs sous une forme elliptique.

En poursuivant l’analyse faite au paragraphe r7.01.32-parametres-elastique-lineaire , juste à l’apparition de l’endommagement, la contrainte longitudinale dans le béton vaut:

\(\frac{{E}_{b}}{1-{\nu}_{b}^{2}}({\epsilon}_{xx}(1-{\nu}_{b}{\nu}_{m})+{\kappa}_{xx}\frac{h}{2}(1-{\nu}_{b}{\nu}_{f}))\)

En confrontant ce résultat à la limite de fissuration du béton \({\sigma}_{t}^{b}\) , on constate que l’on obtient un domaine élastique dans le quadrant \(({N}_{xx},{M}_{xx})\) positifs de forme polygonale:

\(\frac{{N}_{xx}}{{N}_{D}}+\frac{{M}_{xx}}{{M}_{D}}=1\)

On sait que cet endommagement naissant est suivi immédiatement dans un modèle 3D de la plaque en béton armé par l’apparition d’une zone endommagée sur une bonne part de l’épaisseur.

Cette différence de prédiction avec le modèle GLRC_DM est inhérente au choix opéré en (4897) d’une énergie de plaque élasto-endommageable en membrane-flexion, combiné au seuil (3682).

Comme cette différence peut varier entre 0% et 30%, on suggère d’abaisser les valeurs numériques de \({N}_{D}\) et \({M}_{D}\) définies précédemment de 10%.

Bilan de l’identification des paramètres du modèle GLRC_DM

On fait le bilan des expressions analytiques simplifiées proposées dans les paragraphes ci-dessus, exploitant des données géométriques et matériau du béton armé (voir leurs définitions aux sections r7.01.32-parametre-traction à r7.01.32-cas-traction-flexion) servant à établir les valeurs à donner aux mots-clés du modèle GLRC_DM dans code_aster .

Paramètre	mot-clé	Identification	Expression analytique proposée
\({E}_{\mathrm{eq}}^{m}\)	E	module de Young effectif en membrane (unités: force/surface)	\({E}_{\mathrm{eq}}^{m}={E}_{a}\frac{{S}_{a}}{h}+{E}_{b}\cdot \frac{{E}_{b}h+{E}_{a}{S}_{a}}{{E}_{b}h+{E}_{a}{S}_{a}(1-{\nu}_{b}^{2})}\)
\({\nu}_{m}\)	NU	coefficient de Poisson effectif en membrane	\({\nu}_{m}={\nu}_{b}\cdot \frac{{E}_{b}h}{{E}_{b}h+{E}_{a}{S}_{a}(1-{\nu}_{b}^{2})}\)
\({E}_{\mathrm{eq}}^{f}\)	EF	module de Young effectif en flexion (unités: force/surface)	\({E}_{\mathrm{eq}}^{f}=\frac{3}{h}{E}_{a}{S}_{a}{\chi }_{a}^{2}+{E}_{b}\cdot \frac{{E}_{b}h+3{E}_{a}{S}_{a}{\chi }_{a}^{2}}{{E}_{b}h+3{E}_{a}{S}_{a}{\chi }_{a}^{2}(1-{\nu}_{b}^{2})}\)
\({\nu}_{f}\)	NUF	coefficient de Poisson effectif en flexion	\({\nu}_{f}={\nu}_{b}\cdot \frac{{E}_{b}h}{{E}_{b}h+3{E}_{a}{S}_{a}{\chi }_{a}^{2}(1-{\nu}_{b}^{2})}\)
Remarque: ces valeurs peuvent être amendées pour privilégier les situations de chargement de cisaillement plan pur, voir (4904) et (4287).
\({N}_{D}\)	NYT	seuil en traction pure à l’apparition de l’endommagement (unités: force/longueur)	\({N}_{D}={\sigma}_{b}^{t}\frac{{E}_{\mathrm{eq}}^{m}h}{{E}_{b}}\cdot \frac{1-{\nu}_{b}^{2}}{1-{\nu}_{b}{\nu}_{m}}\)
\({\gamma}_{\mathrm{mt}}\le 1\)	GAMMA_T	paramètre de dégradation de raideur en traction	voir r7.01.32-parametre-raideur-traction
\({N}_{C}\)	NYC	seuil en compression pure à l’apparition de l’endommagement (unités: force/longueur)	voir r7.01.32-parametres-compression
\({\gamma}_{\mathit{mc}}\le 1\)	GAMMA_C	paramètre de dégradation de raideur en compression	voir r7.01.32-parametres-compression
\({\alpha}_{c}>0\)	ALPHA_C	Paramètre de retard d’apparition de l’endommagement en compression	voir r7.01.32-parametres-compression
\({M}_{D}\)	MYF	seuil en flexion pure à l’apparition de l’endommagement (unités: force) dans le cas UTIL	voir r7.01.32-parametres-flexion
\({\gamma}_{f}\le 1\)	GAMMA_F	paramètre de dégradation de raideur en flexion	voir r7.01.32-parametres-flexion
Remarque: les valeurs de \({N}_{D}\) et \({M}_{D}\) peuvent être réduites pour limiter l’écart sur la frontière du domaine élastique pour les chargements mixtes en traction-flexion, voir r7.01.32-cas-traction-flexion . On peut aussi vérifier l’évaluation de \({N}_{D}\) sur une situation de cisaillement plan pur, cf. r7.01.32-cas-distorsion .

On pourra avantageusement utiliser l’opérateur DEFI_GLRC, voir [bib12] pour obtenir l’identification des paramètres du modèle GLRC_DM à partir des données des matériaux, acier et béton, et de la géométrie de la section de béton armé.

Il est utile de confronter ces estimations à la réponse donnée par un autre modèle de comportement –comme le modèle ENDO_ISOT_BETON– sur un cas simple, voire d’opérer un recalage sur des courbes de réponse, dans l’intervalle \(\left[0,\mid {\varepsilon}_{xx}^{\max}\mid \right]\) estimé dans l’étude en vue, par exemple en se basant sur le cas-test de vérification SSNS106 [bib8].

Vérification

Ce modèle est vérifié par les tests SSNS106A, B, C, D, E, F (voir [bib8]), par comparaison avec une modélisation multi-couche exploitant le comportement ENDO_ISOT_BETON et des nappes d’acier élastiques. Les cas étudiés sont:

ssns106 a	2 trajets de chargement traction-compression puis compression-traction \(\mid {\varepsilon}_{xx}^{\max}\mid =0,0002\)
ssns106 b	2 trajets de chargement à la suite (flexion -+ puis flexion +-), double cycle
ssns106 c	trajet de chargement combiné avec cyclage en traction 2 fois plus rapide qu’en flexion
sns106 d	cycle de cisaillement pur
ssns106 e	cycle de cisaillement et flexion combinés
ssns106 f et g	cycles de traction-compression et de cisaillement pur avec le kit_ddi GLRC_DM+ VMIS_ISOT_LINE
ssns106 h	traction – compression pure, sollicitations élevées \(\mid {\varepsilon}_{xx}^{\max}\mid =0,001\)
ssns106 i	flexion pure alternée, sollicitations élevées
ssns106 j	couplage traction/compression et flexion, sollicitations élevées
ssns106 k	compression – traction avec ALPHA_C=100
ssns106 l	cisaillement pur et distorsion dans le plan, sollicitations élevées
ssns106 m	couplage flexion et cisaillement dans le plan, sollicitations élevées
ssns106 n	flexion anticlastique, sollicitations élevées
ssns106 o	sollicitations chargement thermique
ssns106 p	traction simple
ssns106 q	compression simple

Validation

On pourra consulter l’ensemble des cas-tests dits de validation dans le domaine du calcul parasismique en transitoire non-linéaire de structures et de bâtiments en béton armé répertoriés dans [bib11].

Bibliographie

[bib1]

LEMAITRE J., CHABOCHE J.L.: «Mécanique des matériaux solides», Ed. Dunod (1985)

[bib2] P.Koechlin, S.Potapov. «A global constitutive model for reinforced concrete plates». ASCE J. Eng. Mech. 2006.

[bib3] P.Koechlin, S.Moulin. «Modèle de comportement global des plaques en béton armé sous chargement dynamique en flexion: loi GLRC améliorée: modélisation de la fissuration par endommagement». Note HT-62/02/021/A, 11/2002.

[bib4] F.VOLDOIRE. «Homogénéisation des structures hétérogènes». Note EDF/DER/MMN HI-74/93/055, 27/10/1993.

[bib5]

S.Moulin. «Modélisation des structures en béton armé sous chargement sismique». Note HT-62/04/025/A, 12/2004.

[bib6]

S.Moulin. F. VOLDOIRE «Étude d’une poutre en béton armé sous chargement de flexion». Note HT-62/05/013/A, 9/2006.

[bib7]

J-J.MARIGO. «Étude numérique de l’endommagement». EDF, Bulletin des études et recherches, série C, n°2, pp.27-48, 1982.

[bib8] (1,2)

[v6.05.106] SSNS106 – Endommagement d’une plaque plane sous sollicitations variées avec la loi de comportement GLRC_DM.

[bib9]

[r3.07.03] – Éléments de plaque DKT, DST, DKQ, DSQ et Q4g.

[bib10] [r7.01.04] – Loi de comportement ENDO_ISOT_BETON.

[bib11]

[a4.01.04] – Dossier de validation du calcul parasismique.

[bib12]

[u4.42.06] – Opérateur DEFI_GLRC.

Annexe : Valeurs propres du tenseur des déformations

On considère une base orthonormée \({({e}_{i})}_{i=1,2,3}\) de l’espace euclidien tridimensionnel, et un tenseur \(\varepsilon\) d’ordre2, symétrique, donc diagonalisable. On note \({\varepsilon}_{j}^{i}\) les composantes mixtes du tenseur \(\varepsilon ={\varepsilon}_{j}^{i}.{e}_{i}\otimes {e}^{\text{*}j}\) dans la base \({({e}_{i})}_{i=1,2,3}\) . On utilise la convention d’Einstein sur les indices mixtes répétés.

On note: \({Q}_{k},{\eta}_{k}\) le triplet des vecteurs propres normés et valeurs propres associées du problème:

(3722)\[\varepsilon .{Q}_{k}={\eta}_{k}{Q}_{k}\iff {\varepsilon}_{j}^{i}.{({Q}_{k})}^{j}{e}_{i}={\eta}_{k}{({Q}_{k})}^{i}{e}_{i}\text{pour}k=1,2,3\]

Remarque

On constate que: \((\varepsilon +\xi \text{Id}).{Q}_{k}=({\lambda}_{k}+\xi ).{Q}_{k},\forall \xi \in ℝ\) , donc ajouter à \(\varepsilon\) un quelconque tenseur diagonal ne modifie pas les directions propres de \(\varepsilon\) .

On sait que les vecteurs propres \({({Q}_{k})}^{j}{e}_{j}\) forment une base orthonormée (repère principal):

(3723)\[{({Q}_{k})}_{i}{e}^{\text{*}i}.{({Q}_{l})}^{j}{e}_{j}={\delta}_{kl}.{\delta}^{{i}_{j}}\Rightarrow {({Q}_{k})}_{j}{({Q}_{\ell})}^{j}={\delta}_{kl}\]

Différencions ces deux relations:

(3724)\[{(d\varepsilon )}_{j}^{i}.{({Q}_{k})}^{j}+{\varepsilon}_{j}^{i}.{(d{Q}_{k})}^{j}=d{\eta}_{k}{({Q}_{k})}^{i}+{\eta}_{k}{(d{Q}_{k})}^{i}\text{pour}k=1,2,3\]

(3725)\[{(d{Q}_{k})}_{j}.{({Q}_{l})}^{j}+{({Q}_{k})}_{j}.{(d{Q}_{l})}^{j}=0\]

Projetons l’équation (3724) sur le vecteur propre \({({Q}_{l})}^{i}{e}_{i}\) et utilisons l’équation (3725):

(3726)\[\begin{split}\begin{array}{}{(d\varepsilon )}_{j}^{i}.{({Q}_{k})}^{j}{({Q}_{l})}_{i}+{\varepsilon}_{j}^{i}.{(d{Q}_{k})}^{j}{({Q}_{l})}_{i}=d{\eta}_{k}{({Q}_{k})}^{i}{({Q}_{l})}_{i}+{\eta}_{k}{(d{Q}_{k})}^{i}{({Q}_{l})}_{i}\text{pour}k,l=1,2,3\\ \iff {(d\varepsilon )}_{j}^{i}.{({Q}_{k})}^{j}{({Q}_{l})}_{i}+{\eta}_{l}.{(d{Q}_{k})}^{j}{({Q}_{l})}_{j}=d{\eta}_{k}.{\delta}_{kl}+{\eta}_{k}{(d{Q}_{k})}^{i}{({Q}_{l})}_{i}\text{pour}k,l=1,2,3\\ \iff {(d\varepsilon )}_{j}^{i}.{({Q}_{k})}^{j}{({Q}_{l})}_{i}=d{\eta}_{k}.{\delta}_{kl}+({\eta}_{k}-{\eta}_{l}).{(d{Q}_{k})}^{i}{({Q}_{l})}_{i}\text{pour}k,l=1,2,3\end{array}\end{split}\]

D’où:

(3727)\[\begin{split}\left\{ \begin{array}{}d{\eta}_{k}={(d\varepsilon )}_{j}^{i}.{({Q}_{k})}^{j}{({Q}_{k})}_{i}\text{pour}k=1,2,3\\ ({\eta}_{k}-{\eta}_{l}).{(d{Q}_{k})}^{i}{({Q}_{l})}_{i}={(d\varepsilon )}_{j}^{i}.{({Q}_{k})}^{j}{({Q}_{l})}_{i}\text{pour}k\ne l=1,2,3\end{array} \right.\end{split}\]

Notons \({\tilde{\varepsilon}}_{j}^{i}\) les composantes mixtes d’un tenseur dans la base \({({Q}_{k})}_{k=1,2,3}\) . Alors:

(3728)\[\begin{split}\left\{ \begin{array}{}d{\eta}_{k}={(d\tilde{\varepsilon})}_{k}^{k}\text{pour}k=1,2,3(\text{pas de sommation sur}k)\\ ({\eta}_{k}-{\eta}_{l}).{(d{Q}_{k})}^{i}{({Q}_{l})}_{i}={(d\tilde{\varepsilon})}_{l}^{k}\text{pour}k\ne l=1,2,3\end{array} \right.\end{split}\]

On vérifie évidemment sur la trace

du tenseur des déformations (qui est indépendante du repère choisi):

(3729)\[\sum_{k=1,2,3}d{\eta}_{k}=\mathrm{tr}(d\tilde{\varepsilon})=\mathrm{tr}(d\varepsilon )=d(\mathrm{tr}\varepsilon )\]

Considérons la densité d’énergie libre d’élasticité isotrope:

(3730)\[\phi (\varepsilon )=\frac{1}{2}\lambda {(\mathrm{tr}\varepsilon )}^{2}+\mu \sum_{k=1,2,3}{({\eta}_{k})}^{2}\]

alors la loi d’état donne le tenseur des contraintes:

(3731)\[\sigma ={\phi }_{,\epsilon }(\varepsilon )=\lambda (\mathrm{tr}\varepsilon )\frac{d\mathrm{tr}\varepsilon }{d\varepsilon }+\mu \sum_{k=1,2,3}{\eta}_{k}\frac{d{\eta}_{k}}{d\varepsilon }=\lambda (\mathrm{tr}\varepsilon )\mathrm{Id}+2\mu \sum_{k=1,2,3}{\eta}_{k}.{({Q}_{k})}^{j}{({Q}_{k})}_{i}.{e}_{j}\otimes {e}^{\text{*}i}\]

En appliquant la remarque faite plus haut, le repère propre du tenseur des contraintes \(\sigma\) est donc identique à celui des déformations \(\varepsilon\) .

Les contraintes principales sont donc naturellement dans le repère principal \({({Q}_{k})}^{j}{e}_{j}\) :

(3732)\[{s}_{k}={\tilde{\sigma}}_{k}^{k}=\lambda \text{tr}(\varepsilon )+2\mu {\eta}_{k}\]

Description des versions du document

Version Code_Aster	Auteur(s) Organisme(s)	Description des modifications
8.4	D.Markovic EDF-R&D/AMA	Texte initial
9.5	S.Fayolle EDF-R&D/AMA	Réécriture des équations et re-formulations de certaines phrases
9.6	F.Voldoire, S.Fayolle EDF-R&D/AMA	Corrections d’équations et re-formulation partielle du modèle; re-rédaction du paragraphe 3 (identification des paramètres). Rédaction de l’annexe: démonstration de la dérivation des valeurs propres.
10.1	F.Voldoire, S.Fayolle EDF-R&D/AMA	Quelques petites corrections et compléments.
10.2	F.Voldoire EDF-R&D/AMA	Modifications de la définition de la position des armatures; ajout de 3 variables internes, voir r7.01.32-variables-internes-modele.
11.1	F.Voldoire, S.Fayolle EDF-R&D/AMA	Ajout des méthodes d’identification pour DEFI_GLRC et du coefficient \({\alpha}_{c}\) .
11.3	F.Voldoire EDF-R&D/AMA	Ajout d’un passage p. 11 expliquant la pente en charge uniaxiale.
13.1	F.Voldoire EDF-R&D/AMA	Compléments pour le cas \({\alpha}_{c}\ne 1\) dans les sections du paragraphe r7.01.32-parametres-elastique-endommageable . Compléments en thermomécanique.
14,3	A. Guilloux EDF-R&D/ERMES	Modifications des paramétrages en flexion et en compression