Homo généis ation numériquede la résistance mécaniquedes massifs fracturés

Une méthode pour déterminer les propriétés derésistance des massifs rocheux fracturés est proposée, àpartir de l'homogénéisation numérique par la méthodedes éléments finis. Les différentes familles de fuacturessont générées à partir des données disponibles et des loisstatistiques. Le comportement mécanique de chaquefracture est modélisé en élastoplasticité, er:r introduisantIes caractéristiques mesurées au laboratoire ou fn situ :

les raideurs normale et tangentielle, la cohésion, l'anglede frottement interne et la dilatance. Le comportement dela matrice est également modélisé en élastoplasticité.4près la détermination numérique du VolumeElémentaire Représentatif, plusieurs simulations sontefFectuées avec différents trajets de charge. Leursrésultats révèlent qu'au-delà du VER, le comportementglobal du massif fracturé est élastoplastique avecécrouissage. Il est caractérisé par une suiface de chargeorientée dont l'expression est donnée. Les valeurs desrésistances anisotropes sont également déterminées.

Mots-clés : massifs rocheux, milieux fracturés,comportement mécanique, critère de résistance,élastoplasticité, homogénéisation.

A. POUYA

M. GHOREYCHIG. 3 S, E c ole p olytechniqyu e

91128 Palaiseau Cedex

_ povya.@q?s.

polytechnique.frrn g @ g 3 s . p o tyt e chnti qu e . fr

l.\r la

I E "/,r..É

| ='rfrI v, ''/#,

l\l Ttn

lÉ'"frr,rj,.|fr.,/li

?;#',/tr)

,fa

T/ti,l'+

/É

,4r,

,r.r,4,

,//ltij

ïrt,ili.t )a4l

.f,.4

Vtl,r.

/./4'/l"t'

|r,/r.r,fr

rtl.,.t/.n

rr,ri,i;/,,,t

14

"/it',.,,,,,1,

rrt'r,ffi.

r,Ii!ïl.tiir.l;

It,IL't(oILl+,I lrll-ot<

NDtË : Les dl'scussrons surcet arttc/e son t acceptéesjusqu'ar.t 'L"' octobre 2001.

Determination of rock massstrength propertiesby numeric al homo geni zation

A method for determining strength properties of fractured rockmasses is proposed. It is based on numerical homogeneisationusing finite element method. Different families of fractures aregenerated on the basis of data available and statistical laws.Mechanical behaviour of each fracture is modelled inelastoplasticity by introducing properties measured inlaboratory or in sifu : oonrâl and tangential stiffnesses,cohesion, angle of internal friction and dilatancy. Matrixbehaviour is also modelled in elastoplasticity. Afterdetermination of Representative Elementary Volume, severalsimulations are performed with different loading paths. Theresults show that beyond REV, global behaviour of fracturedrock mass is elasoplastic with hardening. It is characterised byan oriented criterion the expression of which is given. Thevalues of anisotropic strength are also determined.

Key words; rock mass, fracture, strength criterion, mechanicalbehaviour, elastoplasticity, homogenization.

49REVUE FRANçAIsE oe cÉorucHNteuE

N. 941e'trimestre 2001

EIntroduction

On considère un massif fracturé dont les propriétésmécaniques peuvent être supposées homogènes à uneéchelle suffisamment grande. A cette échelle, le com-portement du massif peut être représenté par celuid'un milieu continu homogène équivalent (MHE). Plu-sieurs méthodes ont été proposées dans la littératurepour la cc classification des massifs fracturés >. Ellespermettant de déterminer, d'une manière empirique etplutôt qualitative, les propriétés homogénéisées dumassif, en fonction de celles de la matrice rocheuse etde certaines caractéristiques de la fracturation (Bartonet aI., 1974 ; Bieniawski, 1976). Lorsqu'une ou deuxdirections préférentielles de la facturation sont consi-dérées, et ceci engendre une périodicité dans la struc-ture du milieu, certaines propriétés mécaniques dumassif peuvent être calculées d'une manière rigoureuseet par des méthodes semi-analytiques, en fonction depropriétés de la matrice rocheuse et des fractures(Bekaert et Maghous, 1996; de Buhan et Maghous,1997).

La démarche présentée ici se situe entre ces deuxextrêmeS : nous considérons le cas de massifs dont lafracturation est exprimée par des lois statistiques quel-conques décrivant l'orientation, l'extension, la densitéet l'épaisseur des fractures. Nous tentons de détermi-ner les propriétés du MHE par la voie numérique. Nousétudions pour cela le comportement d'un domaine suf-fisamment grand du massif fracturé pouvant représen-ter le comportement moyen de ce massif. Le comporte-ment de ce volume élémentaire représentatif (VER) estétudié par la méthode des éléments finis, eh assimilantles fractures à des éléments n joints >. Cette méthode a

déjà été appliquée à l'étude des propriétés d'élasticiténon linéaire du massif (Coste,1997; Coste et a1.,1999).Ce travail en constitue une extension et porte essentiel-lement sur la détermination de la résistance mécaniquedu milieu fracturé dont le comportement global est dutype élastoplastique.

-

Présentation de la méthode

WDescription de la fracturation

Le relevé géologique de la fracturation fournit desinformations sur l'orientation (pendage et azimut), lalongueur, l'espacement (densité de fractures par mètresde carotte) et I'épaisseur des fractures. Ces informationssont souvent très incomplètes. Elles permettent cepen-dant, moyennant certaines hypothèses, de définir desfamilles de fractures qui peuvent être assimilées à desdisques, dans un espace à trois dimensions. Cefle hypo-thèse est fréquemment adoptée dans la modélisation del'écoulement en milieu fracturé (Long et a1.,1987) et peutêtre également admise pour une modélisation méca-nique. Toutefois, nous allons nous restreindre ici à unemodélisation plus simplifiée à deux dimensions. Noussommes conscients de l'insuffisance de cette représen-tation dans Ie cas de milieux fracturés et plaçons cettepremière étude dans la perspective plus large d'une

extension ultérieure à trois dimensions. Ce choix estmotivé par l'objectif spécifique de cette étude quiconsiste à déterminer la résistance mécanique du milieufracturé. Cela exige le recours à de nombreux calculs enélastoplasticité, en présence de nombreuses fractures.Ce type de modélisation n'est pas encore tout à fait opé-rationnel à trois dimensions.

Aussi, dans cette étude, nous allons considérer queles différentes familles de disques sont recoupées pardifférents plans, ce qui permet d'en déduire desmodèles de fractures à deux dimensions.

Dans l'exemple du massif étudié, les fractures ontune extension variant de 20 cm à 30 m suivant une loide distribution statistique exponentielle. Leurs orienta-tions varient entre 0 et 90" par rapport à I'axe x du VER(Fig. 1) et suivent une loi normale avec une moyenned'environ 60" (fractures sub-verticales) et un écart-typede 10o. La densité surfacique des fractures (nombre decentres de fractures par unité de surface) suit une dis-tribution uniforme, avec une valeur moyenne d'environ0,4 m-2 L'épaisseur des fractures est supposée del'ordre de 1 mm (fractures essentieliement colmatées).

ffiComportement de la matrice etdes fractures

Pour simplifler le modèle, nous considérons que lemassif a un comportement mécanique élastoplastiquedu type Mohr-Coulomb. Ce type de comportement estfréquent en mécanique des roches et a servi en partl-culier aux autres approches citées ci-dessus.

On suppose par ailleurs que ]e comportement élas-toplastique de la matrice rocheuse (sans fracture) estcaractérisé par une élasticité linéaire et isotrope et uneplasticité parfaite. On néglige ainsi la rupture fragile dela matrice qui n'intervient que pour des contraintes trèsé1evées par rapport à la résistance au cisaillement desfractures. Le comportement des fractures assimilées àdes cc joints > est également linéaire dans la phase élas-tique et obéit au critère de Mohr-Coulomb dans laphase plastique. Ces comportements simplifiés per-mettent de comparer facilement le comportementhomogénéisé avec celui des constituants, et d'endéduire éventueliement des résultats plus généraux.

Les valeurs des caractéristiques considérées pour lamatrice et les fractures sont les suivantes :

. Roche cristalline :

Élasticifé: module d'Young E = 72 000 MPa, coefficientdePoissonv-0,25Plasticité : cohésion Cinterne ô - 57"

. Fractures naturelles colmatées avec de la calcite :

Elasticité : ratdeur normale k" - 4.105 MP a/m, raideurtangentielle k, : 3.106 MPa/mPlasticité : cohésion c - 1,5 MPa,, angle de frottementinterne e = 27"

WMéthode numérique

Considérons le cas d'un solide avec des fractures(Fig. 1). On impose à la frontière de ce solide des dépla-

50REVUE FRANçAISE DE GEOTECHNIQUEN'941er trimestre 2001

E

Volume ElémentaireReprésentatif

Division du \IER en

sous domaines fermés

Principe de la méthode numérique développée : division du VER en sous-domaines fermés.Numerical method: subdivision of the Representative Elementary Volume (REV) in closed subdomains.

X

cements U' correspondant à des déformations macro-scopiques E,, homogènes, c'est-à-dire de la forme U, -E,,X,, où x, est la coordonnée dans un repère cartésién,où bien des forces Fi correspondant à des contraintesmacroscopiques I,, homogènes, c'est-à-dire de la formeF, = I,,r,, où n est Iâ normale extérieure à la frontière dudomaine au point consid éré. On suppose que ce char-gement crée dans le solide à l'échelle locale (celle deshétérogénéités), des contraintes et déformations dési-gnées par o et e. On note respectivemenl lhom et Eho* lamoyenne (volumique) de o et de e dans le domaine.Notre objectif est d'établir un ( modèle de comporte-ment homogénéisé > reliapl lhom ) phom. Dans le cas oùdes déplacementS U, (resp. forces F,) correspondant àune déformation homogène E (resp. contrainte homo-gène I) sont imposés à la frontière du modèle, on peutmontrer que phom - E (resp. Ihom = I).

Les quantités moyennes lhom ) fhom sont calculées surun volume suffisamment grand (VER) pour être à l'abrides fluctuations statistiques et représentatives du MHE.La taille du VER est déterminée par des calculs méca-niques préliminaires, à partir de la méthode des élémentsfinis, sur des domaines de tailles différentes : le VER cor-respond à une cerfaine taille critique au-delà de laquellela valeur d'un paramètre de référence, par exemple lemodule de Young, re varie plus avec la taille du domaine.Pour l'exemple étudié ici, la taille du VER est d'environ10 m. Le VER correspond donc à un carré (distributionuniforme des fractures) de 10 m de côté (Fig. 1).

Dans Ia pratique, le calcul aux éléments finis néces-site une discrétisation du domaine à partir d'unmaillage respectant la géométrie des fractures. Lesfractures (les joints) ne doivent pas traverser les élé-ments finis (la matrice) du maillage. Pour satisfaire cetteexigence, nous avons dû développer un logiciel per-mettant de diviser automatiquement le VER en sous-domaines fermés. Les fractures en constituent les fron-tières (Fig . 1). Le domaine complet formé par lessous-domaines fermés ainsi créés peut alors être dis-crétisé à l'aide d'un préprocesseur du maillage habi-

tuel. La dernière étape précédant le calcul consiste àdédoubler les næuds correspondant aux fractures,pour former les éléments a joints l.

Le calcul est effectué en contraintes planes ; cettehypothèse facilite la détermination du comportementhomogénéisé dans le plan consid éré

EChargements étudiéset résultats obtenus

Trois types de chargements ont été imposés auVER:1) Compression dans la direction x avec différentesvaleurs de la contrainte latérale ouu

'

9" iryp,oté sur les côtés parallèles à y, oyy constante surles côtés parallèles à x ;

2) Compression dans la direction y avec différentesvaleurs de la contrainte latéral€ o"^ :

[Ju imposé sur les côtés paral]èles à x, oxx constante surle's côtés parallèles à y ;

3) Cisaillement XY avec différentes valeurs de lacontrainte moyenne :

U* imposé sur les côtés parallèles à x, U" imposé surles côtés parallèles à y, (rxx constante sur les'côtés paral-Ièles à x, o* = oxx constante sur les côtés parallèles à y.

Dans chaque cas, les contraintes et 1es déformationshomogénéisées sont calculées à partir de valeurs desforces ou des déplacements imposées ou obtenues surla frontière. Les courbes de contrainte-déformationissues de ces trois types de chargement sont présen-tées sur la figure 2. Pour simplifier l'écriture, sontnotées dans la suite o et r respectivement lescontraintes et les déformations homogénéisées. C'estl'analyse de ces résultats qui permet de déterminer lemodèle du comportement du MHE.

51REVUE FRANçAIsE oe cÉorrcHNteuE

N. g4

l",trimestre 2001

Oxx

Mpa30

25

20

15

10

5

0

0.000 0.001 0.002 0.003 0.004 txx

-rrê^r4 -r4

fr-7 >47-

w' --O- oyy = 0

--I-oyy = lMPa

--+ oyy = 2MPa

rT

Analyse des rësultatsLes résultats ci-dessus peuvent être assimilés aux

cc données expérimentales > qui seraient obtenues lorsd'une expérience de laboratoire, eo vue de l'élaborationd'un modèle rhéologique. On se retrouve confronté auxproblèmes habituels de Ia détermination d'un modèleà partir de données disponibles : celles-ci ne suffisentpas, à elles seules, à Ia mise au point du modèle et doi-vent être complétées par des hypothèses sur la formede la loi de comportement et sur la nature des critèresde résistance. Cela peut apparaître surprenant dans lamesure où la simulation numérique permet a priorid'imposer tout trajet de chargement dans toute direc-tion souhaitée, ce qui signifie qu'il est en principe pos-sible de déterminer rigoureusement le modèle de com-portement sans aucune hypothèse supplémentaire. I1

n'en demeure pas moins qu'en dépit du progrès impor-

i#i#Ëffiiliiit Courbes d'écrouissage du MHErelatives aux différentes directionsde chargement.Stress-strain curves for a homogenizedfractured rock mass corresponding todifferent load directions.

tant des moyens de calcul les simulations numériquesnécessaires à cette étude sont encore longues et fasti-dieuses (plusieurs jours de calcul sur une station SUNSolaris pour une seule simulation) et peuvent difficile-ment être répétées en grand nombre. C'est pourquoinous allons nous contenter de quelques hypothèsescomplémentaires permettant de pousser plus loinl'interprétation des résultats obtenus. Le bien-fondédes hypothèses adoptées pourrait faire l'objet d'unevérification ultérieure.

Cette démarche conduit à proposer un modèlehomogénéisé caractérisé volontairement par desexpressions simples qui impliquent des approximationssur l'évolution des résultats obtenus. En revanche, cer-tains aspects du comportement sont déterminés rigou-reusement par la voie analytique, à partir des proprié-tés des différents constituants. C'est Ie cas par exemplede la limite d'élasticité initiale du MHE, comme nous Ieverrons dans la suite.

5?REVUE FRANçA|sE or cÉorectNteuEN'94ler trimestre 9001

ovv

Mpa

50

40

30

20

10

0

0.000 0.001 0.002 0.003 0.004 Evv

)

L

^6{ -{- oxx = 0

--f- oxx = 10MPa

-+- oxx = 20MPa

f.{pf

0xy

Mpa

I

6

4

2

0

0.0000 0.0008 0.0016

--a-- 0m = 0--F om = lMPa+orn= 2MPa

EModélisation rhéologiq ue

Trois phases de comportement distinctes sontvisibles sur les courbes cc contrainte-déformation l pré-sentées ci-dessus : l'élasticité, la plasticité avec écrouis-sage et la plasticité parfaite. Il convient de soulignerque le MHE présente un écrouissage bien que le com-portement des constituants du milieu soit parfaitementplastique et ne soit donc caractérisé par aucun écrouis-sage. Ce fait traduit un effet de structure du milieuhétérogène. Ce phénomène est également rencontrépour d'autres types de milieux hétérogènes tels que lespolycristaux (Pouya et a1.,1996).

Dans la suite, seront traitées d'abord les deuxphases élastique et parfaitement plastique qui sont plussimples à modéliser . La phase avec écrouissage qui faitla transition entre les deux autres sera étudiée à la fin.

WPhase élastique

Les courbes a contrainte-déformation r commen-cent par une phase élastique Iinéaire s'étendant jusqu'àdes déformations de l'ordre de 10-5 (la Iimite d'élasticitédu massif sera analysée plus loin).

Le tableau I présente les rapports constants entreles différentes variables caractérisant cette phase à par-tir des courbes de compression simple et de cisaille-ment pur.

iitiliiili;:+i*iii|* i*itii Rapports des différentes variables,obtenus dans la phase élastique linéaire.

on remarque que tout chargement dans la directionXX ou suivant YY entraîne des déformations r.,, trèsfaibles, voire théoriquement nulles si le matériau hbmo-généisé est isotrope ou si les directions X et Y sont lesdirections principales de l'anisotropie. Les valeurs der*uÆ** et e*u/e,,u restent d'ailleurs très faibles, négligeablesdëvant celles" du rati o E.^,/E-., définissant le coefficient dePoisson. En outre, les #alâurs des modules d'élasticitéobtenues dans les deux directions XX et YY ne diffèrentpas significativement. Le comportement élastique peutdonc être considéré, en première approximation,comme isotrope. Ce résultat va en faveur d'une hypo-thèse du comportement élastique isotrope et permetd'accéder aux valeurs des caractéristiques élastiqueshomogénéisées du massif, données ci-deSSous :

Ces résultats conduisent à une valeur du module decisaillement G - E/2(1, + v) - 24,7 GPa. C'est précisé-ment la valeur obtenue à partir d'une sollicitation encisaillement (troisième ligne du tableau ci-dessus, oùI'on trouve G - 6,,,/28",, = 24,7 Gpa). Cette vérificationconfirme la cohérëncé''des résultats et la précision descalculs.

ffiLimite de la phase élastique

La détermination de la limite du domaine d'élasti-cité initial est importante, car elle indique le début del'endommagement du massif. Cette limite peut être cal-culée rigoureusement à l'aide d'un calcul analytiqueutilisant Ia théorie du calcul à la rupture. Le domained'élasticité du milieu homogénéisé est l'intersection desdomaines d'élasticité de ses constituants. Comme dansIe cas étudié, la cohésion C de la matrice est très grandedevant celle des fractures c (respectivement 1T et1,51MPa), ce sont les fractures qui atteignent, les pre-mières, la plasticité et qui délimitent, comme nousallons le voir, Ie domaine d'élasticité du massif.

Considérons une fracture faisant un angle 0 avec ladirection x. Soient oo, o* et o* les composantes du ten-

seur des contraintes o. î = (cos0, sinO) et n = (- sinO,cosO) sont les vecteurs respectivement tangent et nor-mal à la fracture. La contrainte de cisaillement t dansla direction t et la contrainte normale o" sur la frac-ture sont donnéeS pâr :

O.,= n.O. n =O-*S

1r \om = ;(o"" + ow /a'

r- t .o.n --ssin20+oxvcosZO

1r \t=rlo* -o*I

En nous plaçant dans le plan (s, o*.,), et en notantdans ce plan S = (s, o".,), U = (co szçii sin20) et V =(- sin20, cos20), Ia condition (1) s'écrit :

cos20 - o*usin?}

rrtË*iifi*iii Échelle localed'une fracture.Local scaie for afracture.

Le domaine d'élasticité de Ia fracture est défini parl'inégalité suivante (la contrainte de compression estnégative) :

Module d'Young :

Coefficient de PoisSoo :

E - 62,6 GPa

v - 0,267

lrl .c-ontgq

lS Vl..-omtgq+S.utgq

s3REVUE FRANçAIsE or cÉorrcuNteuE

N. g4

1"'trimestre 2001

Direction de chargement o*"/r*" [GPa) r/eYY XXt*r/E**

XX 63,2 0,269 0,015

Direction de chargement orr/r* (GPa) r/eXX }rys/exy yy

YY 62,0 0,264 0,005

Direction de chargement o^u/r*u (GPa) t"*/t"u r/tyy xy

XY 49,4 - 0,014 - 0,03 nt\;/

-Â.0ft

Cette condition définit un domaine compris entredeux droites D., et D, faisant un angle q avec la direc-tion U (Fig. 4).'Le dômaine élastique du système des

fractures est l'intersection des domaines ainsi obtenuspour différentes valeurs de 0. Lorsque l'angle 0 variedans un intervalle l0r, 0rl, ce domaine est délimité par

des arcs de cercles et des segments de droites. Cedomaine peut être ouvert ou fermé suivant les valeursde 0.,, 0r. Pour le massif étudié,, 0 varie entre 0 et 90. Dece fait,, l'intersection des domaines élémentaires est un

cercle de centre O et de rayon R donné par :

R - loal sin<p

R-ccosç-o,.sinq

ffiwPhase de plasti citê, parfaite -critàre de résistance

Comme le montre la figure 2, la courbe d'écrouis-sage tend vers une asymptote horizontale exprimantune plasticité parfaite. La contrainte asymptotiquereprésente la résistance maximale du matériau dans ladirection de chargement consid érée. Les valeurs decette contrainte en fonction de la contrainte latérale oude la contrainte movenne sont données dans letabieau II.

Valeurs de la résistance dans différentesdirections en fonction des contrainteslatérales ou movennes.

Direction XX

Contrainte latérale o", (MPa) 0 1 2

Résistance oxx (MPa) 26,0 31.1 35,4

Direction YY

Contrainte latérale o,* (MPa) 0 10 15

Résistance oyy (MPa) 1,4,2 49,0 66,5

Direction XY

Contrainte moyenne o_ (MPa) 0 1 10

Résistance oxv (MPa) 5,1 5,55 11,2

Les résultats rassemblés dans Ie tableau II mettenten évidence une forte anisotropie du comportementplastique, liée aux écarts entre les valeurs de la résis-tance à la compression uniaxiale relatives aux direc-tions x et y. C'est un aspect important du comporte-ment mécanique qui doit être exprimé par le modèledu comportement plastique. Les résultats montrentégalement une dépendance linéaire de la résistancemaximale correspondant aux différentes directions enfonction de la contrainte movenne ou de la contraintelatérale (Fig 5).

(Mpa)

60

50

40

30

20

10

0

0 2 4 6 I 10 12 t+1Upa)

o*, fonction de o'"

--r- o* fonction de la contrainte latérale (orr)

-+ o' fonction de la contrainte latérale (o*r)

;til;iii;;;i1iiliiii:ti1i;'i:,i1|;11ii1|li|iliËIii;$l1ul Effet de la contrainte normale sur larésistance en cisaillement et de la contraintelatérale sur la résistance en compressionsimple dans les directions x et y.Effect of lateral or of mean stress on the shear orthe compressive strength in different directions.

soit :

st &'

lriiilii:Ïlililiilii:llliiii:tiillliiiiiiin ti ljiili Domaine d'élasticité initial du massiffracturé dans une vision homogénéisée.Initial elastic domain of the homogenized rockIT)ASS.

Le domaine élastique du massif est f intersection dece domaine avec celui de Ia matrice. Ce dernier corres-pond à un cercle de centre O avec un rayon donné parla même expression (2) mais caractérisée par C et 0 dela matrice ; le rayon correspondant est donc plus grandque R. En définitif, le domaine élastique du massif estcaractérisé par le cercle de rayon R, tandis que le cri-tère de plasticité initial s'écrit :

o1" +s2 <R'

Cette inégalité peut être exprimée sous la forme :

olu <R=CCoSg

-o"" *ow2 omslnq

On y reconnaît l'expression du critère de Mohr-Coulomb de cohésion c et d'angle de frottementinterne q.

En conclusion, le domaine élastique initial dumassif fracturé homogénéisé est isotrope etdéterminé uniquement par le comportement desfractures.

54REVUE rRANçArsE oe cÉorecuNteuEN" 941er trimestre 2001

Ces résultats permettent de calculer les valeurs del'angle de frottement interne et de la cohésion pour Iemilieu homogénéisé. Le tableau III présente ces valeursobtenues à partir d'une analyse basée sur le charge-ment dans les directions XX, YY et XY.

iiiliiiiiiiiii#iiiilr.*t$ffiiiitiiii:li:: Valeurs de la cohésion et de l'angle defrottement interne homogénéiséscalculées pour les différentes directionsde chargement.

Il convient de remarquer que les valeurs de la cohé-sion et de l'angle de frottement interne homogénéisésse situent bien, comme cela doit être le cas, entre lesvaleurs des fractures (c = 1,51MPa,, A - 27") et les carac-téristiques de la matrice rocheuse (C = 17 }r{rPa, Q - 57").

A partir de ces données, le critère de résistancemaximale correspondant à la phase de plasticité par-faite du MHE a été modélisé moyennant quelqueshypothèses sur la forme de ce domaine et des approxi-mations sur les valeurs des paramètres.

Puisque les valeurs de l'angle de frottement trou-vées dans les directions X et Y sont très proches, ellesseront désormais considérées comme identiques etcaractérisées par la valeur moyenne de 40". Enrevanche, Ies valeurs de la cohésion dans les deuxdirections demeureront bien distinctes.

Supposons maintenant que o,,,, = 0 et plaçons-nousdans le plan (o-, s) (Fig. 6). Les deùx valeurs de Ia cohé-sion données dans le tableau III et la valeur de l'anglede + 40" définissent deux droites dans ce plan délimi-tant le domaine de résistance du milieu (noter quechanger X en Y revient à changer s en - s et donc ç et cen - q et - c). Ces droites passent par les points B = (0,6,02) et C(0, - 3,35) et se croisent au point A(5,59 ; 1,335).Ainsi est obtenue la coupe du domaine de résistancedans le plan (o_, s).

Le domaine de résistance dans le plan (o,,,,s) est délimité par deux droites de pente'l40' par rapport à l'ax€ - o_.The strength domain in (o'_,, s) plane is boundedby two straight lines making an angle + 40" withon., axis.

Fixons maintenant la valeur de o,o et plaçons nousdans le plan (s, o'-) passant par cefte vâieur de o. (Fig. z).Nous pouvons alôrs identifier dans le plan (s, o".,), troispoints appartenant à la frontière du domaine dé résis-tance du milieu : deux points dont les coordonnées sontdonnées par o*u = 0 et par les valeurs de s exprimées parles deux droiteé ci-dessus et un point s = 0 avec o..., = S,0Z- o.tg (31,7").Nous faisons alors 1'hypothèse''que ledomaine de résistance est une ellipse dans ce plan pas-sant par ces points, d'axes parallèles aux axes s et o.^. etcentiée sur ia même ordonnée s que le point A oë tafigure 6. A f intérieur de ce domaine, on retrouve ledomaine d'élasticité initial de MF{E (cercle centré surl'origine) pour la valeur de o. consid érée.

;liiiiiili1iii;;;;l;iiiiiili1;|iii.;iiiiiiiir.*gÏ;iriii;ilr Le domaine de résistance dans le plan (o,.,,s) est une ellipse contenant le domairi'ed'élasticité initial (cercle centré surl'origine).The strength domain in (o*u, s) plane is an ellipsecontaining the initial elasti"c domain.

La forme elliptique supposée pour le domaine derésistance dans le plan (o*.,, s) et la variation linéaire desaxes de cette ellipse aùec la contrainte moyenneconduisent à rechercher le critère de résistance duMF{E sous la forme :

- d)t - (-o,. * nrO,

où a, b, d et g sont quatre constantes à déterminer.

En posant o",, = 0, on trouve :

(- o,' + g)t - o, Cê qui permet defigure 6 et les valeurs du tableau1,59 MPa, g - 5,59 MPa.

f(o-, s, 0) = (a s - d)2 -déterminer, d'après la(3) :â:cotO(40"1,d=

En posant s = 0, on trouve : f(o,'.', 0, o...,) = (b o*, -d)2 -(-o. + g)t - 0, ce qui permet d'obtênir ''Ë = cotg"(3 1,To).

Cette expression conduit à deux droites exprimant lavariation de oxv en fonction de o_ et s'écrivant respecti-vement :

^ o._=: o,n tg(31,7") + 4,43 MPa et - o*u = - om tg(31,7")

"J

- 2,47 MPa.

La valeur de la cohésion c intervenant dans la pre-mière expression est de 4,43 MPa au lieu de 5,02, valeurdonnée dans le tableau III. C'est moyennant cetteapproximation que le critère de résistance est exprimésous une forme simple donnée par (3). Ce modèle pré-voit, par ailleurs, une limite de cisaillement de 2,4T MPadans Ia direction opposée à celle étudiée. Ce point est àvérifier par des calculs complémentaires.

55REVUE FRANçAISE DE GÉOTECHNIQUE

N'941"'trimestre 2001

riiiiiiii:illiiiÏiiiiiii:ii:iiliiiiiiiiiiiiiiiinffi iiô'iii:i:iii

WWModèle d'écrouissage

Le domaine d'élasticité initial du MHE, cercle centrésur l'origine et de rayon R exprimé par (2), peut semettre également sous la forme (3). Les valeurs desparamètres correspondants seront alors les suivantes :

â=b - 1/sin(27') -2,20 d=0- 2,96 MPa.

g - 1,51 cotg(27')

Nous pouvons alors supposer que le domained'élasticité du MHE est défini systématiquement par uncritère de la forme (3) dont les paramètres a, b, d, gdépendent de l'état du système. Nous supposons plusprécisément que ces paramètres dépendent d'une seulevariable d'écrouissage, notée 1,déterminée par ladéformation plastique du milieu à partir de la relation

E = lle ll u,r.c une valeur initiale nulle. Notons alors a0,- il tl

bo, do, go et à* ,b_ , d_ , g_ respectivement les valeursinitiales et asymptotiques de à, b, d et g correspondantà la frontière élastique initiale et au critère de résistancemaximale. Supposons en plus que la dépendance deces paramètres en fonction de € s'exprime sous laforme :

â = âo + (a_ - ao) 0 - e Bç) (4)

avec des formes analogues pour b,, d,, et g.

Le paramètre B est une constante à déterminer parl'ajustement des courbes d'écrouissage. Ce traitementconduit aux valeurs comprises entre 547 et 2 447. Poursimplifier, une valeur unique de F = 1 000 sera retenuepour tous les paramètres a, b, d et g. Ce choix estmotivé par les résultats présentés sur la figure B. Ils'agit de la courbe d'écrouissage dans la direction y.L'ajustement de cette courbe conduit à une valeur de F

caractérisée par cette valeur et celle obtenue avec unevaleur de pnombre de paramètres du modèie, ce qui revient àconsidérer que l'écart observé sur la figure B est accep-table, compte tenu des incertitudes pesant sur laconnaissance des caractéristiques d'un milieu fracturé.

ovv

Mpa

12

10

I6

4

2

0

-./-,l./t/ I

Numérique

F = 2447

F = 1000

I

o

-

{/

0.006 0.008 toyy

En résumé, le comportement plastique du MHE estdécrit comme suit :

Critère de plasticité :

f(o, €) -.{ (o*, s,, oxy, €) = (a s - d)'+ (b o*u d)t - (-)o','+g)t-d2<0Loi d'écrouissage :

a - B(u" - au) . " €

Ë-p(b"-b,,) ."€à-B(d_-d).-"€

s - F(g- - s,) t-" €

*=ll'll

0.002 0.001

Comparaison entre l'ajustement individuel(B = 2 447) de la courbe d'écrouissage dans ladirection y et l'ajustement avec une valeurdeFunique(p=1000).Comparison between the numerical plot(corresponding to compression in y direction),the individual fitting of this plot, leadrng toB - 2,447, and fitting of this plot with a fixed valueof B - 1,000.

Etatinitial: [- 0; a.= ao; b = bo; d = do; g = go

Valeurs numériques : F = 1 000

âo = 2,20; bo = 2,20; do = 0 ;

â- = 1,1,9; b- = 1,,62; d-= 1,59 Mpa ; 9*

EConclusions et perspectives

Une méthode numérique a été présentée pour esti-mer la résistance mécanique de massifs cristallins frac-turés. Elle s'appuie sur le relevé de I'état de fractura-tion du massif et les caractéristiques mécaniques de lamatrice rocheuse et des fractures. Son originalité tient àl'approche mise en æuvre qui consiste à assimiler lemassif à grande échelle à un milieu homogène équiva-lent. Les caractéristiques de la déformabilité et de Iarésistance de ce milieu ont fait l'objet d'une étude fon-dée sur l'homogénéisation numérique. En mettant àprofit les données disponibles sur un massif cristallin,une méthodologie a été mise au point. EIle a permisd'établir un modèle rhéologique susceptible d'êtreemployé ultérieurement pour des études concrètes. Laméthode proposée ne prétend pas être exempte detoute hypothèse, simplification ou approximation. Nouspensons cependant que, par rapport aux méthodesempiriques et qualitatives utilisées en mécanique desroches pour les massifs fracturés, la méthode proposéeintègre davantage et plus rigoureusement l'ensembledes informations pouvant être rassemblées, pour cetype de milieu, à partir d'investigations in situ et aulaboratoire. Ces données concernent en particulierl'état et les caractéristiques mécaniques des fracturesqui jouent un rôle prépondérant dans la résistance dumassif rocheux.

Les résultats obtenus démontrent l'intérêt et la fai-sabilité de la méthode proposée qui peut être étendue àtrois dimensions, sans trop de difficulté. Il sembledésormais possible d'établir des abaques à partir derésultats des simulations numériques qui pourraientêtre réalisées sur quelques configurations typiques deI'état de fracturation des massifs rocheux. ljne foisvalidé par l'expérience de terrain, ce genre d'outil seraittrès utile à l'étude des massifs fracturés. I1 constitueraitun complément indispensable aux outils de modélisa-tion locale employés, par ailleurs, pour l'étude des phé-nomènes à petite échelle.

9o: 2,96 MPa

5ôREVUE FRANçAISE DE GEOTECHNIQUEN'941e' trimestre 200i

WBarton N., Lien R., Lunde J. (1,974) - r< Engi-

neering classification of rock masses forthe design of tunnel support r. RockMechanics, vol. 6, no 4, p.189-236.

Bekaert A., Maghous S. (1996) - < Three-dimensional yield strength properties ofjointed rock mass as a homogenizedmedium >>. Mechanics of Cohesive Fric-tional Materials, vol. 1, p.1-24.

Bieniawski Z.T. (1976) - < Rock mass clas-sifications in rock engineering >. Procee-dings of the symposium on Explorationfor Rock Engineering, Johannesburg.

Buhan (de) P., Maghous S. (1997) - < Com-

portement élastique non linéaire macro-scopique d'un matériau comportant unréseau de joints r. C.R. Acad. Sci. Paris,t.324, Série II b, p.209-218.

Coste F. (1997) - a Comportement thermo-hydro-mécanique des massifs rocheuxfracturés >. Thèse de doctorat, Ecolenationale des ponts et chaussées.

Coste F., Ghoreychl M., Didry O. (1999) -a Modélisation du comportement hydro-mécanique des massifs rocheux fractu-rés par homogénéisation >. 9' Congrèsmondial de mécanique des roches, Bal-kema,'p. 875-880.

Long J.C.S., Billaux D. (1987) - < From fielddata to fracture network modellins : âûexample incorporating spatial structure t.Water Ressources Res earch, vol. 23, no 7 ,

p.1201-1216.Pouya M., Zaoui A., Nabouli M. (1996) - ( A

micro-macro model for polycrystailinehalite >>. The mechanical behavior of salt.Proc 4tt1 Conf., Ghoreychi et al. eds.Trans Tech Pub, p.129-141.

Su K., Ghoreychi M., Chanchole S. (2000) -< Experimental study of damage in gra-nite >. Geotechnique, vol. L, no 3, p. 235-241.

57REVUE FRANçAIsE or eÉorucuNteuE

N" 94'1"'trrmestre 2001