Modélisation du béton avec ANSYS
Modèles plastiques
Dans la simulation numérique non linéaire, les modèles de matériau développés dans le cadre de la théorie de plasticité sont souvent plus stables que ceux développés dans le cadre de la théorie de l’endommagement. C’est pourquoi, on étudie tout d’abord les modèles plastiques du code ANSYS.
Les modèles « Extended Drucker-Prager » (EDP)
Le modèle de Drucker-Prager et les modèles dérivés du critère de Drucker-Prager sont beaucoup utilisés dans la modélisation des géo-matériaux. Dans ANSYS, les modèles «EDP» résultent d’une extension du modèle de Drucker Prager classique. Leur surface de charge et leur potentiel plastique peuvent avoir une forme linéaire ou curviligne. En plus, ils tiennent compte d’un écrouissage non linéaire isotrope.
Caractéristiques des modèles EDP du code ANSYS
Modèles
EDP avec surface de charge linéaire et hyperbolique [19] Les deux formes de la surface de charge sont les suivantes : – Forme linéaire = + . − (ˆ ) = 0 m y pl F q α σ σ ε ( 3.1) avec 3 2 q = J ; 2 J est le 2e invariant du tenseur déviatorique s ( ) ( ) ϕ ϕ α 3 sin 6sin − = ; ϕ est l’angle de frottement interne ( ) y pl σ εˆ est la limite plastique de matériau pl εˆ est la déformation plastique équivalente – Forme hyperbolique (ˆ ) 0 2 2 F = a + q +ασ m −σ y ε pl = avec a > 0 coefficient lié à la forme de surface ( 3.2) L’écoulement plastique peut être associé ou non associé. En cas de plasticité non associée, le 76 potentiel plastique est respectivement: Q = q +α f .σ m −σ y (εˆ pl) = 0 dans le cas linéaire, et ( 3.3) . (ˆ ) 0 2 2 Q = a + q +α f σ m −σ y ε pl = dans le cas hyperbolique, avec α f ≠ α ( 3.4) Ces modèles peuvent reproduire la plasticité parfaite ou l’écrouissage multi linéaire exprimé par la relation entre ( ) y pl σ εˆ et pl εˆ . b) Modèle EDP CAP [19] Contrairement aux deux modèles précédents, la surface de charge EDP de ce modèle est limitée par deux « caps » : « cap » de traction et « cap » de compression. Les paramètres et les formulations de ce modèle sont exprimés en fonction des invariants du tenseur de contrainte : 1 I , 2 J , 3 J . Figure 3.2 Section méridienne de la surface de charge du modèle « CAP » [19] La section méridienne de la surface de charge mixte est illustrée sur la Figure 3.2. Elle se compose de trois parties : enveloppe de cisaillement, « cap » de compression et « cap » de traction. La formule combinée de la surface de charge est présentée dans l’équation ( 3.5). Ici, YC YT YS Γ, , , sont les fonctions qui caractérisent la forme de la surface de charge. ( ) ( ) , , . ( , , ). ( , ). ( , ) 0 , , ( , , , , ) 1 0 2 2 3 2 1 0 0 1 0 2 0 0 1 2 3 0 0 = Γ − = = ψ σ σ σ σ σ σ J J J Y I K Y I Y I Y K Y I J J K C T S ( 3.5) La fonction de l’enveloppe de cisaillement EDP est présentée dans l’équation ( 3.6). La géométrie de cette surface est illustrée sur la Figure 3.3. Dans le cas où A = 0 , elle exprime 77 exactement le critère de Drucker Prager. Le terme 1 . . I Y Ae β a pour but de courber la surface de Drucker Prager pour la faire conformer à la fois au cas de cisaillement simple et au cas de compression bi axiale: 1 . 1 0 0 1 ( , ) . . Y I S Y Y I I Ae β σ =σ −α − ( 3.6) Figure 3.3 Fonction de l’enveloppe de cisaillement [19] Les fonctions du « cap » de compression et de celui de traction sont données dans les équations ( 3.7), ( 3.8) . Elles sont illustrées sur la Figure 3.4 La plasticité se produit sur le « cap » de compression quand 1 K0 I < et sur le « cap » de traction quand I1 > 0 . ( ) 2 0 0 1 0 1 0 0 0 1 . ( , ) ( , , ) 1 . − = − − σ σ R Y K I K Y I K H K I S Y C C ( 3.7) ( ) 2 0 1 1 0 1 . ,0( ) ( , ) 1 . = − σ σ S Y T T R Y I Y I H I avec H : Heaviside fonction Y RC : taux de K0 − X0 par rapport à ( , ) YS K0 σ 0 Y RT : taux de triaxiale T f par rapport à ,0( ) YS σ 0 ( 3.8) 78 Figure 3.4 Fonctions du « cap » de compression (à gauche) et du « cap » de traction (à droite) [19] La forme de la section transversale est déterminée par l’équation ( 3.9). Cette fonction prend en compte la dépendance de la limite plastique en fonction du déviateur de contrainte et du troisième invariant. Cette dépendance est illustrée clairement sur la Figure 3.5. ( ) ( ) ( ) Γ = + + − β ψ β ψ β 1. sin 3 1 . 1 sin 3 2 1 , avec = − − 2/3 2 1 3 .2 .3 .3 sin 3 1 J J β ψ est le taux de la résistance de traction tri axiale par rapport à la résistance de compression ( 3.9) Figure 3.5 Section transversale de la surface de charge [19] Ce modèle permet aussi de modéliser l’écrouissage de la surface de charge. Cela est fait d’une part, avec la loi ( ) p σ 0 γ qui exprime la variation de la résistance de cisaillement simple σ 0 en fonction de la déformation plastique déviatorique effective p γ et d’autre part, avec la loi 79 qui exprime la variation de la résistance de compression tri axiale X0 en fonction de la déformation plastique volumique p v ε selon la formule ( 3.10) : ( ( ))( ) .( 1) 1 2 0 0 = 1 − − − i − i c c p c D D X X X X v ε W e avec c W1 déformation plastique volumique maximale possible Xi résistance compressive tri axiale initiale c c D1 D2 , sont les constants à déclarer du modèle ( 3.10) 3.1.1.2 Evaluation a) Modèle EDP avec la surface de charge linéaire Avec la forme linéaire, la surface de charge n’a que deux paramètres σ Y et α . Elle ne peut pas bien simuler tous les états de contrainte à la fois. On considère tout d’abord la traction uni axiale et la compression uni axiale. Pour que le modèle simule bien la résistance dans ces cas, il faut : 0 3 − . − Y = c c f f α σ ( 3.11) 0 3 + . − Y = t t f f α σ ( 3.12) Ces équations permettent de trouver : . 3 Y 1 c f α σ = − ( 3.13) c t c t f f f f + − = (3 ) α ( 3.14) Pour le béton, t c f = 10 à 6 % f . Ainsi, la valeur deα varie de 45,2 à 66,2 et la valeur de l’angle de frottementϕ varie de 60° à 70° . Malheureusement, avec ces valeurs, la surface de charge ne peut pas modéliser la résistance en compression bi axiale. En effet, l’état de contrainte dans ce cas est suivant: − − = 0 0 0 0 0 0 0 cb cb σ σ σ ( 3.15) 80 cb ⇒ q = .3 J 2 = σ , 2 3/ σ m = − σ cb La surface de charge dans ce cas est suivante: . 0 3 2 -1 0 3 .2 . − = ⇔ − − = cb Y Y cb cb σ σ α σ σ σ α ( 3.16) Evidemment, avec α = 45,2 66,2 à , la partie gauche de l’équation ( 3.16) n’est jamais égale à zéro. C’est-à-dire, la plasticité en compression 2D ne se produit jamais, ce qui n’est pas vrai en réalité. Cela peut être expliqué de façon géométrique : quand l’angle de frottement ϕ = 60° 70 à ° , la surface de cône de Drucker Prager n’a pas d’intersection avec la plan (σ1 , σ 2 , 0 ) dans l’espace des contraintes principales. La plasticité n’a donc pas lieu quand les contraintes sont en état 2D. Si on veut que le modèle simule bien le comportement 2D du béton en compression, l’équation ( 3.16) nous permet de trouver la condition de α : 0 5,1 3 .2 -1 > ⇔ < α α ( 3.17) La valeur raisonnable de α est souvent de 1,2 et l’angle de frottement correspondant est de 30° . Cependant, ces valeurs entraînent une grande valeur de la résistance de traction. Par exemple, si f 60MPa c = , la résistance de traction f t =17MPa , qui est fausse pour le béton. En conclusion, ce modèle n’est pas convenable pour le béton. b) Modèle EDP avec la surface de charge hyperbolique Contrairement au modèle précédent, la forme hyperbolique de ce modèle dispose d’un paramètre de plus a pour courber le sommet du cône Drucker-Prager. Cela permet d’obtenir une résistance t f petite sans imposer une grande valeur de α . Les conditions de charge en compression et traction uni axiales sont respectivement, . 0 3 2 2 F = a + f c − f c −σ y = α et ( 3.18) . 0 3 2 2 F = a + f t + f t −σ y = α ( 3.19) 81 Figure 3.6 Section méridienne de la surface de charge EDP hyperbolique Les équations ( 3.18) et ( 3.19) permettent de déterminer les deux paramètres σ Y et a en fonction de c f , de t f et de α : .2 .( ) .( ) 3 3 2 2 2 c t c t Y f f f f α − − − = α σ ( 3.20) 2 2 . 3 Y c c a f − f = + α σ ( 3.21) La question restante est de trouver la valeur deα pour que la surface de charge modélise bien d’autres états de contrainte. Considérons la compression bi axiale : la condition de charge est la suivante: . 0 3 2 2 2 F = a + f cb − f cb −σ y = α ( 3.22) Comme discuté dans la partie (a), il faut tout d’abord queα ≤ 5,1 . Et puis, avec chaque valeur deα , les paramètresσ Y et a sont déterminés à l’aide des équations ( 3.20), ( 3.21). La résistance cb f est ensuite calculée grâce à l’équation ( 3.22). Tous ces calculs ont été faits dans Excel pour deux types du béton : le premier a f c = 60MPa et f t = 4MPa , le deuxième a f c = 30MPa et f t = 3MPa . La résistance en compression bi axiale calculée est représentée sur la Figure 3.7. On trouve que la résistance en compression bi axiale modélisée est beaucoup plus grande que la valeur réelle. L’erreur minimale de calcul possible est de 58% 72 114 72 = − pour le 1e cas et de 58% 36 57 36 = − pour le 2e cas.