4. Lois de comportement pour les bétons

La liste suivante concerne les lois de comportement pour le béton.

4.1. Modèle de Mazars

'MAZARS' : Modele d'endommagement scalaire pour le béton (bien adapté aux chargements monotones).

4.1.1. Description

La loi de comportement élastique endommageable en traction/compression selon le modèle de Mazars, corrigé pour prendre en compte de manière plus réaliste l'endommagement en cisaillement, est présentée dans [MAZARS-1984] [MAZARS-1986] [PIJAUDIER-1991].

Caracteristiques et limitations principales :

  • Dans ce modèle la dégradation des propriétés élastiques du matériau est représentée à l'aide d'une variable scalaire \(D\) variant entre zéro (matériau sain) et l'unité (matériau totalement endommagé). Cette dernière est obtenue par la combinaison de deux variables scalaires représentant l'endommagement sous sollicitations de compression et de traction séparément ;

  • Cela permet de modéliser convenablement la dyssimétrie traction-compression observée expérimentalement pour les bétons. Cependant, aucune reprise de raideur induite par la refermeture des fissures lors du passage d'une sollicitation de traction à une sollicitation de compression ne peut être prise en compte, c'est-à-dire que l'effet unilatéral n'est pas modélisé ;

  • Compte tenu de cela, cette loi est adaptée à la simulation de la réponse du béton sous chargement monotone, mais nécessite des modifications pour une utilisation dans le cadre d'un calcul sous sollicitations cycliques et/ou dynamiques. Dans ce dernier cas, des limitations supplémentaires sont présentes (par exemple, l'impossibilité de modéliser des boucles d'hystérésis).

4.1.2. Formulation du modèle

La relation contrainte-déformation s'écrit :

\[\boldsymbol{\sigma} = (1-D) \mathbb{E} : \boldsymbol{\varepsilon}\]

où \(\boldsymbol{\sigma}\) est le tenseur des contraintes de Cauchy, \(\boldsymbol{\varepsilon}\) est le tenseur des déformations infinitésimales, et \(\mathbb{E}\) est le tenseur d'élasticité d'ordre 4.

La variable d'endommagement évolue au cours du chargement en fonction du critère d'endommagement :

\[f = f(e) = e - \kappa\]

où \(e\) est une mésure scalaire du tenseur de déformation et \(\kappa\) est une variable d'histoire fournissant le maximum atteint par la déformation équivalente pendant le chargement du matériau.

Selon la formulation proposée par Mazars, la déformation équivalente est définie par :

\[{e}=\sqrt{\langle\boldsymbol{\varepsilon}\rangle:\langle\boldsymbol{\varepsilon}\rangle} = \sqrt{\sum_{i=1}^{^{n}}\langle\epsilon_{i}\rangle^{2}}\]

où \(\langle\boldsymbol{\varepsilon}\rangle\) est la partie positive du tenseur des déformations, \(\epsilon_{i}\) est la i-ème déformation principale, \(\langle\cdot\rangle_+\) est l'opérateur de Macaulay, et \(n\) répresente la dimension du problème consideré.

La variable d'histoire \(\kappa\) est donc définie par :

\[\kappa = \max_t (\kappa,e) \qquad \kappa(t=0) = e_0\]

avec \(e_0\) une valeur seuil initiale et \(t\) une variable répresentant le temps (ou bien le pseudo-temps dans un calcul quasi-statique).

L’évolution de cette surface de charge doit respecter les conditions Kuhn-Tucker :

\[f \leq 0 \qquad \dot{\kappa} \geq 0 \qquad \dot{\kappa} f = 0\]

En d'autres termes, la surface de charge ne peut croître que si le seuil de déformation actuel est dépassé. Cela implique que l'endommagement ne progresse pas pendant les phases de décharge ou les phases de charge où les niveaux de déformation sont inférieurs au maximum atteint précédemment au cours de l'historique du chargement.

Pour prendre en compte la nature fortement dissymétrique du comportement en traction et en compression du béton, Mazars introduit deux fonctions \(D_t\) et \(D_c\) représentant respectivement les dégradations en traction et compression. Elles sont définies comme suit :

\[D_{t(c)} = 1 - \frac{e_0 (1-A_{t(c)})}{\kappa} - A_{t(c)} \exp\left[ -B_{t(c)}(\kappa - e_0)\right]\]

avec \(A_{t(c)}\) et \(B_{t(c)}\) les quatre paramètres additionnels permettant de définir, avec le seuil de première fissuration en traction \(e_0\), les lois d'évolution de l'endommagement en traction (t) et en compression (c). Le paramètre \(A_{t(c)}\) permet de controler la contrainte résiduelle en traction (respectivement compression) uniaxiale tandis que le paramètre \(B_{t(c)}\) contrôle la forme de la loi d'evolution de l'endommagement dans la phase post pic de contrainte.

La variable d'endommagement \(D\) est finalement obtenue par combinaison linéaire des variables \(D_{t}\) et \(D_{c}\) comme suit :

\[D = \alpha_t^\beta D_t + \alpha_c^\beta D_c\]

avec \(\alpha_{t(c)} \in [0,1]\) des facteurs de combinaison qui s'expriment en fonction des déformations principales comme suit :

\[\alpha_t = \frac{\sum_{i=1}^{n} \varepsilon_i^t \langle \varepsilon_i \rangle_+}{e^2} \qquad \alpha_c = 1 - \alpha_t\]

avec \(\varepsilon_i^t\) les déformations associées aux contraintes principales positives. Le paramètre \(\beta\) a été introduit historiquement plus tard dans le modèle pour éviter une évolution trop rapide de l'endommagement en cisaillement [PIJAUDIER-1991].

4.1.2.1. Réponses typiques

../_images/Figure_Mazars_1.png

Critère de Mazars. (a) Surface seuil dans l'espace des contraintes. (b) Trace dans le plan \(\sigma_3=0\). Giry 2001.

../_images/Figure_Mazars_2.png

Loi contrainte - déformation pour une sollicitation uniaxiale.

4.1.2.2. Quelques commentaires

Grâce à sa simplicité et sa robustesse, ce modèle a été et est encore largement utilisé pour modéliser le comportement du béton. Certaines pathologies peuvent néanmoins être citées et pour lesquelles des développements sont à considérer :

  • Le modèle présente une fragilité excessive dans son comportement en cisaillement et l'introduction du paramètre \(\beta\) pour atténuer cet effet entraîne une reprise de rigidité à des niveaux de déformation élevés ;

  • Le modèle ne prend pas en compte l'effet unilatéral, c'est-à-dire une reprise de raideur due à la refermeture des fissures expérimentalement observée. En conséquence, le modèle ne parvient pas à reproduire correctement le comportement sous chargements cycliques ;

  • En termes numériques, l'utilisation de l'opérateur de Macaulay dans l'expression des coefficients \(\alpha_{t(c)}\) entraîne une dérivée non définie de ceux-ci en zéro. Cela empêche ainsi l'utilisation de l'opérateur tangent dans le schéma de résolution. Par conséquent, seul l'opérateur sécant est utilisé, ce qui limite la vitesse de convergence du schéma de résolution ;

  • Le caractère isotrope de l’endommagement ne permet pas de bien suivre l’évolution des nonlinéarités pour des chargements non radiaux.

4.1.3. Implémentation Cast3M (esope)

Dans la suite, nous détaillons les étapes du calcul pour les éléments volumiques (sources cmazar.eso) d'une part et pour les éléments poutres à fibres (sources fibmaz.eso) d'autre part, en mettant l'accent sur les parties de code correspondantes aux aspects théoriques mentionnés précédemment. Pour une analyse détaillée de l'implémentation et des aspects plus strictement techniques concernant la signification des variables, veuillez vous référer aux commentaires présents dans les fichiers sources.

4.1.3.1. Implémentation pour les éléments volumiques

L'implémentation est réalisée dans le fichier source cmazar.eso

1C CMAZAR    SOURCE    MB234859  25/03/18    21:15:01     12209          
2      SUBROUTINE CMAZAR (WRK52,WRK53,WRK54,WRKK2,NSTRS1,NVARI,
3     1                   ICARA,JDIM,IFOUR2)
4.1.3.1.1. Entrées
10C     WRK0     pointeur sur un segment deformation au pas precedent
11C
12C     WRK1     pointeur sur un segment increment de deformation
13C
14C     WRKK2     pointeur sur un segment variables internes au pas precedent
15C
16C     WRK5      pointeur sur un segment de deformations inelastiques
17C
18C     XMATER     constantes du materiau
19C
20C     NSTRS1      nombre de composantes dans les vecteurs des contraintes
21C                                        et les vecteurs des deformations
22C
23C     NVARI      nombre de variables internes (doit etre egal a 2)
24C
25C     NMATT      nombre de constantes du materiau
26C
27C     ISTEP      flag utilise pour separer les etapes dans un calcul non local
28C                ISTEP=0 -----> calcul local
29C                ISTEP=1 -----> calcul non local etape 1 on calcule les seuils
30C                ISTEP=2 -----> calcul non local etape 2 on continue le calcul
31C                               a partir des seuils moyennes

35C     JDIM       Dimension de travail
36C                ( Coques JDIM =2 , Massifs JDIM = IDIM )
37C     IFOUR2     Type de formulation
38C                ( Coques IFOUR2 = -2 => contraintes planes ,
39C                  Massifs IFOUR2 = IFOUR)
4.1.3.1.2. Sorties
43C     VARINF      variables internes finales
44C
45C     SIGMAF      contraintes finales
4.1.3.1.3. Algorithme

Le calul de l'endommagement est réalisé par une procédure purement explicite.

  • On calcule la déformation totale au niveau du point d'intégration ;

  • On calcule le tenseur des déformations principales ;

  • On calcule la matrice d'élasticite et les contraintes principales ;

  • On calcule la déformation équivalente de Mazars :

    • Si le calcul est local (ISTEP = 0), la déformation principale est évaluée directement sur la base des déformations principales ;

    • En cas d'un calcul non-local, l'évolution de l'endommagement est pilotée par la contrepartie non-locale de la déformation de Mazars. Celle-ci est évaluée avec deux passages dans la loi de comportement :

      • Lors du premier passage (ISTEP = 1), on calcule la déformation locale et on sort de la loi de comportement. La déformation non-locale est calculée via une procédure ad-hoc en dehors de la loi de comportement, par exemple, via une méthode non-locale intégrale ou bien une formulation de type gradient implicite ;

      • Cette déformation non-locale est une variable d'entrée de la loi de comportement (ISTEP = 2) et est utilisée pour faire évoluer l'endommagement ;

  • On vérifie le dépassement du seuil de déformation. Si le seuil n'est pas dépassé, l'endommagement n'est pas mis à jour. Sinon, on procède comme suit.

  • On calcule les coefficients \(\alpha_{t(c)} \in [0,1]\). Pour cela :

    • On calcule le signe des contraines elastiques de compression SIGPC(i) (négative) et de traction SIGPT(i) (positive) et les traces associées :

      184            DO 300 ISTRS=1,3
185               IF (SIGP(ISTRS).LT. XZERO) THEN
186                  SIGPC(ISTRS)=SIGP(ISTRS)
187                  SIGPT(ISTRS)=XZERO
188               ELSE
189                  SIGPT(ISTRS)=SIGP(ISTRS)
190                  SIGPC(ISTRS)=XZERO
191               END IF
192300         CONTINUE
193            TRSIGT=SIGPT(1)+SIGPT(2)+SIGPT(3)
194            TRSIGC=SIGPC(1)+SIGPC(2)+SIGPC(3)

    • On calcule les déformations associées aux contraintes positives \(\varepsilon_i^t\) :

      198            DO 400 ISTRS=1,3
199               EPSILT(ISTRS)=(SIGPT(ISTRS)*(UN+XNU)-TRSIGT*XNU)/YOUN
200400         CONTINUE

    • On calcule \(\alpha_{t}\) puis on en déduit \(\alpha_{c}\) :

      204            ALFAT = MAX(XZERO,EPSIPP(1))*EPSILT(1) +
205     1              MAX(XZERO,EPSIPP(2))*EPSILT(2) +
206     2              MAX(XZERO,EPSIPP(3))*EPSILT(3)
207            ALFAT = ALFAT/(EPSTIL*EPSTIL)
208            ALFAC = UN - ALFAT

    • On corrige les paramètres de combinaison linéaire via le coefficient \(\beta > 1\) pour amémiorer la réponse en cisaillement :

      221            IF (BETA .GT. UN) THEN
222               IF ( ALFAT .GT. XPETIT ) THEN
223                  ALFAT=ALFAT**BETA
224               END IF
225               IF ( ALFAC .GT. XPETIT ) THEN
226                  ALFAC=ALFAC**BETA
227               END IF
228            END IF

  • On corrige la déformation equivalente de Mazars \(e\) par le coefficient \(\gamma\) pour améliorer la réponse en bi ou tri-compression :

    \[e = e \gamma\]

    le coefficient \(\gamma\) est calculé de la façon suivante :

    \[\gamma = -\frac{\sqrt{\sum_{i=1}^n \langle \sigma_i \rangle_{-}^2}}{\sum_{i=1}^n \langle \sigma_i \rangle_{-}}\]

    avec \(\langle \cdot \rangle_{-}\) l'opérateur partie négative.

    212            IF (TRSIGC.LT. -YPETIT .AND. TRSIGT.LT.YPETIT) THEN
213               GAMMA=SIGPC(1)*SIGPC(1)+SIGPC(2)*SIGPC(2)+
214     1                                 SIGPC(3)*SIGPC(3)
215               GAMMA=-SQRT(GAMMA)/TRSIGC
216               EPSTIM=EPSTIM*GAMMA
217            END IF

  • Le calcul de la variable d'endommagement D est effectué après avoir vérifié si le seuil de dommage initial a été dépassé. Cette vérification est nécessaire car il est possible que la valeur de la déformation equivalente de Mazars ait été multipliée par \(\gamma\). Tandis que l'évolution du dommage en compression DC suit le modèle de Mazars classique, trois lois d'évolution différentes du dommage en traction DT sont proposées selon la valeur du paramètre ATRA :

    • ATRA > 0 : modèle de mazars classique ;

    • -10 < ATRA < 0 : loi d'évolution exponentielle modifiée pour prendre en compte l'énergie de fissuration \(G_{f}\) via le paramètre BTRA ;

    • ATRA < -10 : loi d'évolution linéaire. Le paramètre BTRA représente alors la déformation pour laquelle la contrainte s'annule.

    252            IF (EPSTIM .GT. EPSD0) THEN
253               IF (ATRA .GT. 0.D0) THEN
254                 DT=UN - EPSD0*(UN-ATRA)/EPSTIM -
255     1              ATRA*EXP(BTRA*(EPSD0-EPSTIM))
256               ELSE IF (ATRA . GT. -10.) THEN
257                 DT=UN - epsd0/epstim*EXP(BTRA*(EPSD0-EPSTIM))
258               ELSE
259                 IF(EPSTIM .LT. BTRA) THEN
260                   DT=UN - EPSD0*(BTRA - EPSTIM)/EPSTIM/(BTRA - EPSD0)
261                 ELSE
262                   DT=1.D0
263                 ENDIF
264               END IF
265               DC=UN - EPSD0*(UN-ACOM)/EPSTIM -
266     1            ACOM*EXP(-BCOM*(EPSTIM-EPSD0))
267            ELSE
268
269               DT=XZERO
270               DC=XZERO
271            END IF
272            D = ALFAT*DT + ALFAC*DC

    La variable d'endommagement est ensuite bornée supérieurement à 0.99999999 afin d'éviter un trop mauvais conditionnement de la matrice de rigidité ;

  • On calcule la nouvelle contrainte et on sort de la loi de comportement ;

  • Les données de sortie sont la contrainte et les variables internes mises à jour.

4.1.3.2. Implémentation pour les éléments poutres à fibres

L'implémentation est réalisée dans le fichier source fibmaz.eso.

1C FIBMAZ    SOURCE    MB234859  24/07/09    21:15:03     11959          
2      SUBROUTINE FIBMAZ (XMAT,DEPS,SIG0,VAR0,SIGF,VARF)
4.1.3.2.1. Entrées
14C     XMAT : Caracteristique du materiaux
15C
16C     EPS : Deformations initiales
17C
18C     DEPS : Increment de deformations
19C
20C     SIG0 : Contraintes initiales
21C
22C     VAR0 : Variables internes initiales
4.1.3.2.2. Sorties
27C     SIGF : Contraintes finales
28C
29C     VARF : Variables internes finales
4.1.3.2.3. Algorithme

Le calul de l'endommagement est réalisé selon la même méthode que celle décrite précédemment dans le fichier source cmazar.eso pour les éléments volumiques, dans la mesure où il s'appuie sur la formulation 3D complète du modèle de Mazars.

  • On calcule la déformation totale au niveau du point d'intégration ;

  • On calcule les terme de la matrice 3x3 de déformation à partir de la théorie des poutres de Timoshenko ;

    65      EPS33(1,1) = VARF(3)
66      EPS33(1,2) = (VARF(4)/2.D0)
67      EPS33(2,1) = EPS33(1,2)
68      EPS33(1,3) = (VARF(5)/2.D0)
69      EPS33(3,1) = EPS33(1,3)
70      EPS33(2,2) = ((-1.D0*XNU)*VARF(3))
71      EPS33(3,3) = EPS33(2,2)
72      EPS33(3,2) = XZERO
73      EPS33(2,3) = XZERO

  • On calcule le tenseur des déformations principales ;

  • On calcule les contraintes principales à partir des déformations principales de de la matrice de Hook 3D (cas général) ;

    86      SIGP(1)=((UN-XNU)*EPSIPP(1)+
87     1         (XNU*EPSIPP(2))+
88     2         (XNU*EPSIPP(3)))
89      SIGP(2)=((UN-XNU)*EPSIPP(2)+
90     1         (XNU*EPSIPP(1))+
91     2         (XNU*EPSIPP(3)))
92      SIGP(3)=((UN-XNU)*EPSIPP(3)+
93     1         (XNU*EPSIPP(2))+
94     2         (XNU*EPSIPP(1)))
95      DO 200 ISTRS=1,3
96         SIGP(ISTRS)=(SIGP(ISTRS)*YOUN/((UN+XNU)*(UN-2.D0*XNU)))
97200   CONTINUE

  • On calcule la déformation équivalente de Mazars en se plaçant toujours dans le cas où le calcul est local (ISTEP = 0). En conséquence, elle est évaluée directement sur la base des déformations principales ;

    102      EPSTIL=MAX( XZERO , EPSIPP(1) )**2 +
103     1       MAX( XZERO , EPSIPP(2) )**2 +
104     2       MAX( XZERO , EPSIPP(3) )**2
105      EPSTIL=SQRT (EPSTIL)
106C
107C      on est toujours dans le cas istep=0
108C
109         EPSTIM=EPSTIL
110         VARF(1)=EPSTIL

  • On calcule les coefficients \(\alpha_{c(t)} \in [0,1]\). Pour cela :

    • On calcule le signe des contraines elastiques de compression SIGPC(i) (négative) et la trace associée :

      120            DO 300 ISTRS=1,3
121                  SIGPC(ISTRS)=MIN(XZERO,SIGP(ISTRS))
122300         CONTINUE
123            TRSIGC=SIGPC(1)+SIGPC(2)+SIGPC(3)

    • On vérifie le dépassement du seuil de déformation. Si le seuil n'est pas dépassé, l'endommagement n'est pas mis à jour. Sinon, on procède comme suit ;

    • On calcule les déformations associées aux contraintes négatives \(\varepsilon_i^c\) :

      142            DO 400 ISTRS=1,3
143               EPSILC(ISTRS)=(SIGPC(ISTRS)*(UN+XNU)-TRSIGC*XNU)/YOUN
144400         CONTINUE

    • On calcule \(\alpha_{c}\) puis on en déduit \(\alpha_{t}\) :

      148            ALFAC = MAX(XZERO,EPSIPP(1))*EPSILC(1) +
149     1              MAX(XZERO,EPSIPP(2))*EPSILC(2) +
150     2              MAX(XZERO,EPSIPP(3))*EPSILC(3)
151            ALFAC = ALFAC/(EPSTIL*EPSTIL)
152            ALFAT = UN - ALFAC

    • On corrige les paramètres de combinaison linéaire via le coefficient \(\beta > 1\) pour amémiorer la réponse en cisaillement :

      156            IF (BETA .GT. UN) THEN
157               IF ( ALFAT .GT. XPETIT ) THEN
158                  ALFAT=ALFAT**BETA
159               END IF
160               IF ( ALFAC .GT. XPETIT ) THEN
161                  ALFAC=ALFAC**BETA
162               END IF
163            END IF

  • A noter que la correction \(\gamma\) en cas de bi ou tri-compression n'a pas de sens avec la modélisation poutre à fibres qui ne traite que des chargements de type traction-compression dans la direction de sa fibre neutre et cisaillement dans le plan de sa section.

  • Le calcul de la variable d'endommagement D est effectué en ayant vérifié au préalable si le seuil de dommage initial a été dépassé. L'évolution du dommage en compression DC comme celle du dommage en traction DT suivent le modèle de Mazars classique :

    169               DT=UN - EPSD0*(UN-ATRA)/EPSTIM -
170     1            ATRA*EXP(BTRA*(EPSD0-EPSTIM))
171               DC=UN - EPSD0*(UN-ACOM)/EPSTIM -
172     1            ACOM*EXP(BCOM*(EPSD0-EPSTIM))
173               D = ALFAT*DT + ALFAC*DC

    La variable d'endommagement est ensuite bornée supérieurement à 0.99999 afin d'éviter un trop mauvais conditionnement de la matrice de rigidité ;

  • On calcule la nouvelle contrainte et on sort de la loi de comportement ;

  • Les données de sortie sont la contrainte et les variables internes mises à jour.

4.1.4. Implémentation MFront

Une implémentation du modèle de Mazars a été réalisée sous MFront par Elian Dussart lors de son stage (2025) dans le cadre du Pôle de Compétences du SEMT. Elle s'inspire du fichier source cmazar.eso utilisé dans Cast3M et reprend la formulation originale de l'endommagement de Mazars sans régularisation mais avec les améliorations de la réponse du modèle en cisaillement et en situation de bi ou tri-compression [DUSSART-2025].

  1//**********************************************************************
  2
  3@DSL Default;
  4@Behaviour Mazars;
  5@Author Elian Dussart;
  6@Date 03 / 06 / 2025;
  7@Description {
  8  Loi de comportement de Mazars.
  9}
 10
 11@UnitSystem SI;
 12
 13//**********************************************************************
 14//********** Material properties to be entered in Cast3M ***************
 15
 16@MaterialProperty stress young;
 17young.setGlossaryName("YoungModulus");
 18@MaterialProperty real nu;
 19nu.setGlossaryName("PoissonRatio");
 20
 21@MaterialProperty real acomp;
 22@MaterialProperty real bcomp;
 23@MaterialProperty real atrac;
 24@MaterialProperty real btrac;
 25@MaterialProperty real beta;
 26@MaterialProperty real epsd0;
 27
 28//**********************************************************************
 29//***************** State (internal) variables *************************
 30//********(updated during calculations and outputted)*******************
 31
 32//@AuxiliaryStateVariable : donnée interne du modèle sauvegardée pendant le pas de temps, permettant d'utiliser la version actualisée dans le pas de temps courant.
 33@AuxiliaryStateVariable real d; // Damage variable (scalar)
 34d.setGlossaryName("Damage");
 35@StateVariable real dcom ;
 36dcom.setEntryName("EndommagementCompression") ;
 37@StateVariable real dtra ;
 38dtra.setEntryName("EndommagementTraction") ;
 39@StateVariable real alft ;
 40alft.setEntryName("ZoneTraction") ;
 41@StateVariable real alfc ;
 42alfc.setEntryName("ZoneCompression") ;
 43
 44//**********************************************************************
 45//********************** Local variables *******************************
 46
 47@LocalVariable real epstim;
 48@LocalVariable real gamma;
 49@LocalVariable real eto_tild;
 50@LocalVariable real Zero;
 51@LocalVariable real Un;
 52@LocalVariable real Deux;
 53@LocalVariable real X_Petit;
 54@LocalVariable real Y_Petit;
 55@LocalVariable real D_Borne;             // maximum endommagement numérique utilisé
 56@LocalVariable real lambda;              // constante de Lamé
 57@LocalVariable real mu;                  // constante de Lamé
 58@LocalVariable stensor<3, stress> sig_p; // contraintes principales
 59@LocalVariable stensor<3, stress> spt;   
 60@LocalVariable stensor<3, stress> spc;   
 61@LocalVariable stensor<3, strain> eto_t; 
 62
 63//**********************************************************************
 64//******** Initialization of Lamé's Constants & Initial Damage *********
 65
 66@InitializeLocalVariables {
 67  lambda = computeLambda(young, nu);
 68  mu = computeMu(young, nu);
 69  D_Borne = 1. - 1.e-8;
 70  Zero = 0.;
 71  Un = 1.;
 72  Deux = 2.;
 73  X_Petit = 1.e-12;
 74  Y_Petit = 1.;
 75}
 76
 77//**********************************************************************
 78//************ Constitutive model integrator (explicit) ****************
 79
 80@Integrator {
 81
 82  const auto e = eval(eto + deto);                                                        // Calcul de la déformation totale
 83
 84  strain e1, e2, e3;       
 85  stress sp1, sp2, sp3;    
 86
 87  constexpr auto esolver = Stensor::FSESJACOBIEIGENSOLVER;                                //Déclaration du solver utilisé
 88
 89  e.template computeEigenValues<esolver>(e1, e2, e3);                                     //Calcul des déformations principales
 90
 91  const auto sig_tilde = eval(lambda * trace(e) * Stensor::Id() + 2 * mu * e);            //récupération des contraintes liées à la déformation calculée.
 92  
 93  sig_tilde.template computeEigenValues<esolver>(sp1, sp2, sp3);                          // Calcul des contraintes principales
 94
 95  sig_p = {sp1, sp2, sp3, stress(0), stress(0), stress(0)};                               // Définition du tenseur de contrainte
 96
 97  eto_tild = sqrt(pow(max(Zero, e1), 2) + pow(max(Zero, e2), 2) + pow(max(Zero, e3), 2)); // Calcul de la déformation équivalente
 98  
 99  eto_tild = max(X_Petit, eto_tild);
100  epstim = max(X_Petit, eto_tild);
101  // Deux variables différentes car la variable epstim peut changer pour la bi ou tri-compression																		  
102  
103  //Pour les calculs structurels, vérification de la condition de déformation (évite des calculs inutiles)
104  if (epstim >= epsd0) {
105	  // Calcul du signe des contraintes élastiques
106	  for (size_t i = 0; i != 3; i++) {
107		if (sig_p(i) < 0) {
108		  spc(i) = sig_p(i);
109		  spt(i) = Zero;
110		} else {
111		  spc(i) = Zero;
112		  spt(i) = sig_p(i);
113		}
114	  }
115	  const auto tr_spc = trace(spc);
116	  const auto tr_spt = trace(spt);
117
118	  // Calcul des deformations dues aux contraintes positives
119	  for (size_t i = 0; i != 3; i++) { // déformation traction
120		eto_t(i) = (spt(i) * (Un + nu) - tr_spt * nu) / young;
121	  }
122
123	  // Déduction alft et alfc
124	  alft = (max(Zero, e1) * eto_t(0) + max(Zero, e2) * eto_t(1) +
125			   max(Zero, e3) * eto_t(2)) /
126			  (pow(eto_tild, 2));
127	  alfc = Un - alft;
128	  
129	  // Modification pour la bi ou tri-compression
130	  if ((tr_spc < -Y_Petit) && (tr_spt < Y_Petit)) {
131		gamma = -sqrt(pow(spc(0), 2) + pow(spc(1), 2) + pow(spc(2), 2)) / tr_spc;
132		epstim *= gamma;
133	  }
134
135	  // Amélioration de la réponse en cisaillement pour beta > 1
136	  if (beta > Un) {
137		if (alft > X_Petit) {
138		  alft = pow(alft, beta);
139		}
140		if (alfc > X_Petit) {
141		  alfc = pow(alfc, beta);
142		}
143	  }
144
145	  // Calcul de dtra et dcom
146	  // Endommagement Mazars classique (at > 0)
147	  if (epstim >= epsd0) {
148		dtra = Un - epsd0 * (Un - atrac) / epstim - atrac * exp(btrac * (epsd0 - epstim));
149		dcom = Un - epsd0 * (Un - acomp) / epstim - acomp * exp(-bcomp * (epstim - epsd0));
150	  } else {
151		dtra = Zero;
152		dcom = Zero;
153	  }
154	  
155	  //Calcul de l'endommagement d
156	  if(alft * dtra + alfc * dcom > d) {
157		d = alft * dtra + alfc * dcom;
158		d = min(d, D_Borne);
159	  }
160  } else {
161	  d = max(d, Zero);
162  }
163
164  // Calcul des contraintes finales
165  sig = (1 - d) * sig_tilde;
166
167}
168
169//**********************************************************************

4.1.5. Hypothèses de calcul et éléments finis disponibles

  • Cette loi est disponible pour les éléments massifs 3D et 2D sous l'hypothèse de contraintes/déformations planes (@lien vers section éléments finis?@).

  • De plus, elle est également applicable aux poutres à fibres (éléments finis poutre à modèle section). Dans ce dernier cas, le modèle a été implémenté dans le modèle poutre à fibres selon sa formulation 3D complète, plutôt qu'uniaxiale.

  • Elle peut être utilisée avec des éléments de type coque sous l'hypothèse de contraintes planes.

4.1.6. Mots clefs dans les opérateurs MODE et MATE

Exemple d'utilisation du modèle de Mazars pour des éléments finis massifs 3D CUB8 :

MODE maillage 'ELASTIQUE' 'ENDOMMAGEMENT' 'MAZARS' 'CUB8' ;
MATE modele 'YOUN' val_youn 'NU' val_nu ('RHO' val_rho)
            'KTR0' val_e0 'ACOM' val_ac 'BCOM' val_bc 'ATRA' val_at 'BTRA' val_bt 'BETA' val_beta ;

Exemple d'utilisation du modèle de Mazars pour des éléments finis poutres à fibres constituée d'éléments finis TIMO dont la section est constituée d'éléments finis QUAS :

  • pour la section :

MODE mail_section 'ELASTIQUE' 'PLASTIQUE' 'MAZARS' 'QUAS' ;
MATE mode_section 'YOUN' val_youn 'NU' val_nu
                  'KTR0' val_e0 'ACOM' val_ac 'BCOM' val_bc 'ATRA' val_at 'BTRA' val_bt 'BETA' val_beta
                  'ALPY' 1. 'ALPZ' 1. ;

  • pour la poutre :

MODE mail_poutre 'ELASTIQUE' 'SECTION' 'PLASTIQUE' 'SECTION' 'TIMO' ;
MATE mode_poutre 'MODS' mode_section 'MATS' mate_section 'VECT' (0. 1. 0.) ;

4.1.7. Paramètres de la loi non linéaire

  • KTR0 : seuil en déformation pour la traction, \(e_0\)

  • ATRA : paramètre pour la traction, \(A_t\)

  • ACOM : paramètre pour la compression, \(A_c\)

  • BTRA : paramètre pour la traction, \(B_t\)

  • BCOM : paramètre pour la compression, \(B_c\)

  • BETA : correction pour le cisaillement, \(\beta\)

4.1.7.1. Valeurs typiques

Pour un béton ordinaire, on peut choisir :

  • \(e_0= 10^{-4}\)

  • \(A_t= 1\)

  • \(A_c= 1,5\)

  • \(B_t= 8000\)

  • \(B_c= 1550\)

  • \(\beta= 1\)

4.1.7.2. Prise en compte de la régularisation dans la définition des paramètres matériaux

  • Régularisation énergetique

  • Régularisation non-locale (quelle formulation? quelle variable est rendue non-locale?)

4.1.8. Anomalies observées et limitation numérique

4.1.8.1. Anomalie 1

Une anomalie a été identifiée dans les sources fibmaz.eso (modèles poutres à fibres) et cmazar.eso (modèles massifs) concernant le mauvais calibrage du test de bicompression dans le modèle de Mazars. Dans fibmaz.eso, le critère trop sévère (1.D-12) générait à tort la correction \(\gamma\) de bicompression, même lorsqu'on est en situation de traction simple ; tandis que, dans cmazar.eso, il conduisait à mal calculer les contraintes en situation de bicompression.

En ce qui concerne la source cmazar.eso pour les modèles massif, ce test à été recalibré à 1 Pa, valeur jugée suffisamment proche de 0 selon le REX (Thèse de Martin Debuisne, 2024).

  • Constatée le 24/07/2024.

  • Anomalie #12209 corrigée dans la version 2025 de Cast3M.

En ce qui concerne la source fibmaz.eso pour les modèles poutres à fibres, la correction \(\gamma\) a été inhibée car la biaxialité du chargement n'a pas de sens avec l'élément poutre à fibres qui ne traite que des chargements de type traction-compression dans la direction de sa fibre neutre et cisaillement dans le plan de sa section.

  • Constatée le 24/07/2024.

  • Anomalie #11959 corrigée dans la version 2025 de Cast3M.

4.1.8.2. Anomalie 2

  • Constatée le 13/06/2024.

  • Anomalie #11948 corrigée dans la version 2024.1 de Cast3M.

Une anomalie a été identifiée dans la source idendo.eso. Elle est datée du 21/08/2023 et cause une erreur d'initialisation du paramètre BETA dans cmazar.eso. Cette anomalie ne rend pas le modèle de Mazars inutilisable mais corrompt ses résultats avec des éléments volumiques. Elle impacte la version 2024.0 de Cast3M.

4.1.8.3. Anomalie 3

  • Constatée le 04/02/2025.

  • Non corrigée à ce jour.

Une anomalie de fonctionnement des modélisations poutre à fibres a été identifiée dans différents cas de chargement élémentaire (traction, compression, flexion, ...) et avec différents modèles de comportement (Mazars, ACIER_UNI, ...). Il s'agit d'une instabilité numérique de type "flambement" qui conduit soit à un arrêt du calcul par non convergence de PASAPAS, soit à des résultats erronés produits après l'instabilité, comme ici un champ de dommage devenant hétérogène dans la section du modèle poutre à fibres, avec une bifurcation entre les résultats de deux paires de points de Gauss conduisant à des valeurs anormales (cf. figure ci-dessous).

../_images/Figure_anomalie3_1.png ../_images/Figure_anomalie3_2.png ../_images/Figure_anomalie3_3.png ../_images/Figure_anomalie3_4.png
../_images/Figure_anomalie3_4.png

Anomalie constatée dans le cas d'une modélisation poutre à fibres soumise à un chargement uniaxial de compression en déplacement imposé à l'extrémité libre d'une poutre encastrée à son autre extrémité - Maillage et conditions aux limites, champ de dommage hétérogène dans la section en fin de calcul, évolutions anormales au cours du chargement de la contrainte moyenne en fonction de la déformation moyenne et du dommage en chaque point de Gauss de la section en fonction du temps.

Cette anomalie, qui concerne donc le modèle poutre à fibres et pas le modèle de Mazars, est contournée dans tous les cas rencontrés jusqu'ici en bloquant les rotations du point correspondant à l'extrémité libre soumise au chargement de la poutre à fibres (cf. figure ci-dessous).

../_images/Figure_solution3_1.png ../_images/Figure_solution3_2.png ../_images/Figure_solution3_3.png ../_images/Figure_solution3_4.png
../_images/Figure_solution3_4.png

Anomalie contournée en bloquant les rotations du point d'extrémité libre de la poutre - Maillage et conditions aux limites, champ de dommage homogène en fin de calcul, évolutions normales au cours du chargement de la contrainte moyenne en fonction de la déformation moyenne et du dommage en chaque point de Gauss de la section en fonction du temps.

4.1.8.4. Limitation numérique

Le modèle de Mazars dans Cast3M, tant dans la configuration éléments volumiques (source cmazar.eso) que poutres à fibres (source fibmaz.eso), exhibe un domaine post-ruine consolidant non physique. Cet artefact numérique est dû à la limitation du dommage maximum :

\[D_{max}=(1 - \epsilon)\]

où \(\epsilon\) est un paramètre arbitrairement petit, défini dans les sources Cast3M du modèle de Mazars, permettant de se prémunir de l'absence complète de rigidité aux points de Gauss ayant atteint la ruine, ce qui empêcherait la poursuite du calcul. L'augmentation de la valeur de ce paramètre est favorable à la stabilité numérique mais défavorable au réalisme de la simulation.

En effet, la consolidation qui en découlerait dans une zone jugée trop grande du modèle E.F. peut conduire à des résultats numériques qui ne sont pas physiquement admissibles et ainsi fausser le jugement du spécialiste du béton, ce qui est préjudiciable à la confiance accordée au modèle.

Le choix de la valeur du paramètre \(\epsilon\) résulte donc d'un compromis à faire entre convergence numérique et représentativité physique du modèle.

L'historique des valeurs attribuées au paramètre \(\epsilon\) est le suivant :

  • dans la configuration éléments volumiques (source cmazar.eso) :

    • \(\epsilon=10^{-20}\) : valeur d'origine ;

    • \(\epsilon=10^{-4}\) : valeur dans la version du jour de Cast3M jusqu'au 9/10/2024 ;

    • \(\epsilon=10^{-8}\) : valeur dans la version du jour de Cast3M depuis le 10/10/2024 ;

  • dans la configuration éléments poutres à fibres (source fibmaz.eso) :

    • \(\epsilon=10^{-5}\) : valeur d'origine.