Lois de comportement pour les bétons

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

Loi MAZARS

'MAZARS' : Modele d'endommagement scalaire pour le beton (bien adapte aux chargements monotones).

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 ici [MAZARS-1984] [MAZARS-1986] [PIJAUDIER-1991].

Characteristiques et limitations principales :

• Dans ce modèlela 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 lors du passage d'une sollicitation de compression à une sollicitation de traction 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).

Anomalies observées

Compte tenu des résultats des cas tests de vérification @, l'utilisation de cette loi dans le cadre d'une modèlisation de type multi-fibre est proscrite !.

Formulation du modèle

La rélation 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 infinitesimales, et \(\mathbb{E}\) est le tenseur d'elasticité d'ordre 4.

La variable d'endommagement evolue au cours du chargement en fonction 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 MacAuley, 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 le 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(\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}\) permet de controler la contrainte résiduelle en traction uniaxiale tandis que le paramètre \(B_{t}\) controle 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 lineaire 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 = \sum_{i=1}^{n} \frac{\varepsilon_i^t \langle \varepsilon_i \rangle_+}{e} \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 tatd dans le modèle pour éviter une évolution trop rapide de l'endommagement en cisaillement [PIJAUDIER-1991].

Réponses typiques

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

Loi contrainte - déformation pour une sollicitation uniaxiale.

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 que la refermeture des fissures expérimentalement observée entraîne une reprise de raideur. 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 Mac Cauley 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.

Implémentation Cast3M (esope)

@Détailler les sources de l'implémentation multi-fibre@

Dans la suite, nous détaillons les étapes du calcul 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 le fichier source mazars.eso.

C MAZARS    SOURCE    CHAT      05/01/13    01:38:20     5004
      SUBROUTINE MAZARS (WRK0,WRK1,WRKK2,WRK5,NSTRS,NVARI,
     1                   NMATT,ISTEP,ICARA,JDIM,IFOUR2)

Entrées

C     WRK0     pointeur sur un segment deformation au pas precedent
C
C     WRK1     pointeur sur un segment increment de deformation
C
C     WRKK2     pointeur sur un segment variables internes au pas precedent
C
C     WRK5      pointeur sur un segment de deformations inelastiques
C
C     XMATER     constantes du materiau
C
C     NSTRS      nombre de composantes dans les vecteurs des contraintes
C                                        et les vecteurs des deformations
C
C     NVARI      nombre de variables internes (doit etre egal a 2)
C
C     NMATT      nombre de constantes du materiau
C
C     ISTEP      flag utilise pour separer les etapes dans un calcul non local
C                ISTEP=0 -----> calcul local
C                ISTEP=1 -----> calcul non local etape 1 on calcule les seuils
C                ISTEP=2 -----> calcul non local etape 2 on continue le calcul
C                               a partir des seuils moyennes

C     JDIM       Dimension de travail
C                ( Coques JDIM =2 , Massifs JDIM = IDIM )
C     IFOUR2     Type de formulation
C                ( Coques IFOUR2 = -2 => contraintes planes ,
C                  Massifs IFOUR2 = IFOUR)

Sorties

C     VARINF      variables internes finales
C
C     SIGMAF      contraintes finales

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 coéfficients \(\alpha_{t(c)} \in [0,1]\). Pour cela faire :

  • On calcule le signe des contraines elastiques :

                DO 300 ISTRS=1,3
                   IF (SIGP(ISTRS).LT. XZERO) THEN
                      SIGPC(ISTRS)=SIGP(ISTRS)
                      SIGPT(ISTRS)=XZERO
                   ELSE
                      SIGPT(ISTRS)=SIGP(ISTRS)
                      SIGPC(ISTRS)=XZERO
                   END IF
    300         CONTINUE
                TRSIGT=SIGPT(1)+SIGPT(2)+SIGPT(3)
                TRSIGC=SIGPC(1)+SIGPC(2)+SIGPC(3)
    

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

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

  • On calcule \(\alpha_{t(c)}\) :

                ALFAT = MAX(XZERO,EPSIPP(1))*EPSILT(1) +
         1              MAX(XZERO,EPSIPP(2))*EPSILT(2) +
         2              MAX(XZERO,EPSIPP(3))*EPSILT(3)
                ALFAT = ALFAT/(EPSTIL*EPSTIL)
                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 :

                IF (BETA .GT. UN) THEN
                   IF ( ALFAT .GT. XPETIT ) THEN
                      ALFAT=ALFAT**BETA
                   END IF
                   IF ( ALFAC .GT. XPETIT ) THEN
                      ALFAC=ALFAC**BETA
                   END IF
                END IF
    

• On corrige la déformation equivalente pour améliorer la réponse en bi- ou tri-compression. Pour cela faire, on modifie \(e\) comme suit :

  \[e = e \gamma \qquad \gamma = \frac{\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.

              IF (TRSIGC.LT. -XPETIT .AND. TRSIGT.LT.XPETIT) THEN
                 GAMMA=SIGPC(1)*SIGPC(1)+SIGPC(2)*SIGPC(2)+
       1                                 SIGPC(3)*SIGPC(3)
                 GAMMA=-SQRT(GAMMA)/TRSIGC
                 EPSTIM=EPSTIM*GAMMA
              END IF
  

• Le calcul de la variable d'endommagement est effectué après avoir vérifié si le seuil initial a été dépassé. Cette vérification est nécessaire car il est possible que la valeur ait été multipliée par \(\gamma\) :

              IF (EPSTIM .GT. EPSD0) THEN
                 DT=UN - EPSD0*(UN-ATRA)/EPSTIM -
       1            ATRA*EXP(-BTRA*(EPSTIM-EPSD0))
                 DC=UN - EPSD0*(UN-ACOM)/EPSTIM -
       1            ACOM*EXP(-BCOM*(EPSTIM-EPSD0))
              ELSE
                 DT=XZERO
                 DC=XZERO
              END IF
              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.

Implémentation MFront

Une implémentation de la loi de Mazars a été réalisée sous MFront. Le code suivant détaille l'implémentation pour une utilisation avec des elements volumiques/surfaciques. La formulation implémentée est une version simplifiée de celle disponible dans Cast3M. En particulier, aucun correctif n'est pas introduit pour améliorer la réponse du modèle en cisaillement et compression bi-/tri-axiale.

//**********************************************************************
@DSL Default;
@Behaviour Mazars;
//@ModellingHypotheses{"PlaneStress"};
@Author Giuseppe Rastiello;
@Date 27 / 10 / 2021;
@Description{Original Mazars model formulation (1984) without successive
			improvements to better account for damage in shear model and
			bitension/bicompression};

//**********************************************************************
//********** Material properties to be entered in Cast3M ***************

@MaterialProperty stress young;
young.setGlossaryName("YoungModulus");  //Young modulus
@MaterialProperty real nu;
nu.setGlossaryName("PoissonRatio");		//Poisson's ratio

@MaterialProperty real Ac;				//Mazars' model parameter
@MaterialProperty real At;				//Mazars' Model parameter
@MaterialProperty real Bc;				//Mazars' Model parameter
@MaterialProperty real Bt;				//Mazars' Model parameter
@MaterialProperty real ed0;				//Damage threshold (strain)

//**********************************************************************
//***************** State (internal) variables *************************
//********(updated during calculations and outputted)*******************

@StateVariable real d;					//Damage variable (scalar)
d.setGlossaryName("Damage");
@StateVariable real kappa;				//History variable
@StateVariable real eeqc;				//Mazars equivalent strain

//**********************************************************************
//********************** Local variables *******************************

@LocalVariable stress lambda;			//Lamé’s constant
@LocalVariable stress mu;				//Lamé’s constant

@InitLocalVariables{
  lambda = computeLambda(young,nu);
  mu     = computeMu(young,nu);
}

//**********************************************************************
//************ Constitutive model integrator (explicit) ****************
@Integrator {

  // total strain tensor (epsilon) and trace of the strain tensor
  const auto e  = eval(eto + deto);
  const auto tr = trace(e);
  
  // eigenvalues of the strain tensor
  strain e1,e2,e3;
  e.template computeEigenValues<Stensor::FSESJACOBIEIGENSOLVER>(e1,e2,e3);
  //without eval at line 50 computeEigenValues does not work !
  
  // Mazars equivalent strain
  const real ppe1=max(strain(0),e1);
  const real ppe2=max(strain(0),e2);
  const real ppe3=max(strain(0),e3);
  eeqc= sqrt(ppe1*ppe1+ppe2*ppe2+ppe3*ppe3);
  
  const auto f = eeqc - kappa ;
  
  if (f > 0  and eeqc > ed0) {
  // Ver 2
  //  if (f > 0) {
  
	//set history variable
	kappa = eeqc ;
	// Ver 2
	//kappa = max(kappa,ed0) ;
	   
	// eigenvalues of the effective stress tensor
	const stress s1 = 2*mu*e1+lambda*tr;
	const stress s2 = 2*mu*e2+lambda*tr;
	const stress s3 = 2*mu*e3+lambda*tr;
  
	// strains due to the positive stresses
	const stress pps1=max(stress(0),s1);
	const stress pps2=max(stress(0),s2);
	const stress pps3=max(stress(0),s3);  
	const stress trpps = pps1+pps2+pps3 ;
	const real elt1 = (pps1*(1.+nu)-trpps*nu)/young;
	const real elt2 = (pps2*(1.+nu)-trpps*nu)/young;
	const real elt3 = (pps3*(1.+nu)-trpps*nu)/young;
  
	// compute alphat and alphac
	const real alphat = (ppe1*elt1+ppe2*elt2+ppe3*elt3)/(eeqc*eeqc) ;
	const real alphac = 1. - alphat ;
    
	// compute dt and dc (note: variable dt cannot be used in Mfront)
	const real dtrac = 1. - ed0*(1-At)/kappa - At*exp(-Bt*(kappa-ed0)) ;
	const real dcomp = 1. - ed0*(1-Ac)/kappa - Ac*exp(-Bc*(kappa-ed0)) ;

	// compute damage
	d = alphat*dtrac + alphac*dcomp;
	d = min(0.99999,d);
		
  }
  // compute stress
  sig = (1.-d)*(lambda*tr*Stensor::Id()+2*mu*e);
}

//**********************************************************************
//************* Elastic/Secant/Tangent stiffness tensor ****************
@TangentOperator{
      if(smt==ELASTIC){
        Dt = lambda*Stensor4::IxI()+2*mu*Stensor4::Id();
      } else if(smt==SECANTOPERATOR){
        Dt = (1.-d)*(lambda*Stensor4::IxI()+2*mu*Stensor4::Id());
      } else {
        return FAILURE;
      }
}      
//**********************************************************************

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 multi-fibres (éléments finis de section). Dans ce dernier cas, le modèle a été implémenté dans le modèle à fibre 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.

Mots clefs dans l'opérateur MODE

Exemple d'utilisation de la loi Mazars pour des éléments finis de section CUB8 :

MODE maillage 'ELASTIQUE' 'ENDOMMAGEMENT' 'MAZARS' 'CUB8' ;

Exemple d'utilisation de la loi Mazars pour des éléments finis de section QUAS :

MODE mail_section 'ELASTIQUE' 'PLASTIQUE' 'MAZARS' 'QUAS' ;

Paramètres de la loi non linéaire

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

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

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

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

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

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

Valeurs typiques

Pour un béton ordinaire, on peut choisir : - \(e_0= 10^{-4}\) - \(A_c= 10^{-4}\) - \(A_t= 10^{-4}\) - \(B_t= 10^{-4}\) - \(B_c= 10^{-4}\)

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

• Régularisation énergetique

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