Généralités sur la modélisation

Au cours de ces dernières années, la modélisation est devenue une des composantes de la méthode scientifique, au même titre que ses deux composantes traditionnelles à savoir la théorie et l’expérience. C’est évidemment l’apparition des ordinateurs puissants qui est à l’origine de ce développement rapide de la simulation numérique. Ce développement a permis de trouver des solutions à de nombreux problèmes, autrefois long, voir irrésolubles, avec les moyens de calcul classiques. En effet il existe aujourd’hui des modèles numériques permettant d’étudier des problèmes aussi complexes que l’action des trafics routiers et ferroviaires, des glissements de terrains, des montées de niveau d’eau, etc. La modélisation est quasiment devenue aujourd’hui un point de passage obligé pour qualifier les structures en configuration extrêmes. Souvent même s’ils sont utilisés, les logiciels de mécanique des sols ne sont pas exploités au maximum de leurs possibilités, ceci à cause de la méconnaissance de leur méthode de fonctionnement et des modèles qu’ils proposent, mais aussi de la difficulté d’interprétation des résultats.

Eléments de mécanique des milieux continus

Par définition une matrice est un tableau rectangulaire de nombres écrit entre crochets :

Les nombres a_ij sont appelés éléments de la matrice et le premier indice (i) d’un élément indique le numéro de la ligne et le second celui de la colonne dans lequel se trouve l’élément.

Une matrice est dite carrée si le nombre de lignes est égal au nombre de colonnes. On note une matrice carrée d’ordre n par (a_nn) et les éléments a₁₁ ;a ₂₂ et a_nn sont appelés éléments diagonaux. La trace d’une matrice A = (a_ij), carrée d’ordre n, est définie comme étant la somme des éléments diagonaux et on note tr(A) :

Si A et B sont deux matrices carrées d’ordre n et l un réel alors :

tr (A + B) = tr A + tr B

tr (AB) = tr (BA)

tr (lA) = l tr A

Deux matrices A= (a_ij) et B= (b _ij) sont dites égales si et seulement si elles sont de même ordre et si chaque éléments de l’une est égal à l’élément correspondant de l’autre c’est à dire a_ij = b_ij. La somme de deux matrices A= (a_ij) et B= (b _ij) est une matrice C= (c_ij) dont les éléments sont définis par c_ij = a_ij + b _ij on remarquera qu’on ne peut additionner deux matrices que si elles sont de même ordre. Si A est une matrice d’ordre (n, p) et B une matrice d’ordre (p, q), on définit la matrice C = AB d’ordre (n, q) par :

La condition d’existence du produit AB est que le nombre de colonne de A doit être égal au nombre de ligne de B. La transposée d’une matrice A d’ordre (n, m) est une matrice B notée ^tA d’ordre (m, n) et de terme général obtenu en échangeant les lignes et les colonnes de A :

"i Î í1,…..,ný et "jÎí1,…..,mý b_ij = a _ji

Propriétés : ^t(^tA)= A ; ^t(lA) = l ^tA ; ^t(A + C) = ^tA + ^tC  ; ^tAC = ^tC ^tA.

Notions de contraintes et de déformations.

Soit un corps (G) en équilibre. Sous l’action des forces extérieures P₁, P₂ , P_n, des forces extérieures prennent naissance et s’exercent entre les différentes parties du corps.

(G)

Forces extérieures s’exerçant sur le corps (G)

Pour étudier la grandeur de ces forces en un point quelconque (O), supposons le corps (G) divisé en deux parties (A) et (B) par une section transversale quelconque (mm) passant par ce point O. En considérant l’une des parties ainsi formées (partie A), on peut admettre qu’elle est en équilibre sous l’action des forces extérieures P₄, P₅, P_n et des forces intérieures réparties sur toute la surface de la section transversale (mm) et qui représente les actions moléculaires exercées par la matière constituant la partie (B) sur celle constituant la partie (A). Considérons par hypothèse que ces forces sont réparties d’une manière continue sur la surface (mm) (répartition homogène). Pour obtenir la grandeur de la quantité de force agissant sur l’élément de surface dA, on admet les actions moléculaires exercées par la partie (B) sur la partie (A) à travers l’air élémentaire dA peuvent être remplacées par leur résultante dP. On appelle contrainte la quantité de force agissant sur l’unité de surface :

Répartition homogène des forces

Lorsque les forces ne sont pas uniformément réparties sur la section transversale, on définit la contrainte en ce point O comme étant la limite du rapport

quand dA tend vers le point O.

Une contrainte a les propriétés d’un vecteur. Sa direction fait en général un angle avec la surface élémentaire sur laquelle elle agit. On peut donc la décomposer en deux composantes :

- la contrainte normale s perpendiculaire à la surface d’application,

- la contrainte tangentielle ou de cisaillement t dans le plan de cette surface.

Représentation des contraintes sur un élément cubique

Les contraintes n’ont de sens que par rapport à une facette que l’on oriente par sa normale. Ainsi si on prend comme élément de volume un cube, les contraintes normales sont les contraintes s_x, s_y et s_z et les indices x, y et z indiquent la direction de la composante elle-même. Les contraintes de cisaillement sont :

t_xy, t_yz, t_zx

t_xz, t_yx, t_zy

Le premier indice indique la direction de la normale au plan de la contrainte et le second celle de la composante elle-même.

L’étude de la répartition des contraintes autour d’un point permet de tirer les conclusions suivantes :

q les extrémités des vecteurs représentant les contraintes décrivent, lorsque l’orientation de la facette sur laquelle s’applique la contrainte varie, un ellipsoïde ayant sont centre au point considéré (ellipsoïde de Lamé) ;

q il existe trois plans orthogonaux privilégiés appelés plans principaux sur lesquels les contraintes sont normales au plan autrement dit l’obliquité ainsi que la composante en cisaillement est nulle. On note en générale ces contraintes en s₁, s₂ et s₃ (avec s₃ £ s₂ £ s₁). s₁ est la contrainte principale majeur, s₂ la contrainte principale moyenne ou intermédiaire et s₃ la contrainte principale mineure.

La notion de tenseur est abstraite. Le tenseur des contraintes est représenté par la matrice carrée :

Il faut donc neuf (9) composantes pour définir le tenseur par contre si le tenseur est symétrique ce qui est le cas, seul six composantes suffisent pour le définir complètement et dans ce cas on a : t_xy = t_yx ; t_xz = t_zx et t_yz = t_zy. Lorsqu’on adopte comme axe de coordonnées les trois directions principales de l’ellipsoïde, le tenseur des contraintes sera défini par la matrice carrée :

Il est toujours possible de décomposer ce tenseur de contraintes principales (diagonales) en un tenseur sphérique et au maximum trois tenseurs déviatoriques :

P est la contrainte moyenne elle est égale à :

S₁ =

S₂ =

S₃ =

Le premier tenseur est un tenseur sphérique qui représente une pression hydrostatique d’intensité P égale à la moyenne arithmétique des trois contraintes principales. Les trois autres tenseurs sont des tenseurs déviatoriques ou déviateur. Remarquons que la trace de chacun de ces déviateurs est nulle et que la somme des trois déviateurs donne le déviateur principal. Le tenseur des contraintes présente trois invariants dont les expressions sont dans l’espace des contraintes quelconque :

1. I₁ = s_x + s_y + s_z

2. I₂ = s_xs_y + s_ys_z +s_zs_x -t²_xy - t²_yz - t²_zx

3. I₃ = s_xs_ys_z + 2t_xyt_yzt_zx -s_xt²_yz -s_yt²_zx - s_zt²_xy

Si on considère comme axe de coordonnée les trois axes principaux autrement si on travaille dans l’espace des contraintes principales, les expressions des invariants du tenseur s’écrivent sous la forme :

1. I₁ = s₁ + s₂ + s₃

2. I₂ = s₁s₂ + s₂s₃ + s₃s₁

3. I₃ = s₁s₂s₃

Notons en fin que le tenseur des contraintes est symétrique c’est à dire que s_ij = s_ji et qu’il représente une approximation assez grossière de ce qui se passe en chaque point du matériau.

Notion de déformation

Définition de la déformation d’un milieu continu

Les déformations sont engendrées par des sollicitations extérieures mécaniques, par des évolutions chimiques et des réorganisations structurelles, etc. Elles posent un problème de compatibilité à savoir elles sont compatibles avec la continuité du milieu. Les déformations imposées, non compatibles, induisent des efforts intérieurs qui, s’ils sont excessifs, peuvent entraîner des désordres, voir même des ruptures. Après la description des contraintes, nous allons introduire ici le tenseur des déformations, notée, plus ou moins complexe en chaque point du matériau. On considère, afin d’étudier le comportement du milieu à l’occasion de la variation du système de charge qu’il supporte et de décrire la déformation, un petit élément de volume entourant un point M. On sait qu’une déformation infinitésimale de ce petit volume est le produit d’une translation, d’une rotation et d’une déformation pure (distorsion).

Soit deux points M et M’appartenant au volume V.

Déformation infinitésimale du milieu

Dans la déformation qui transforme le volume V en un volume V₁, ces points viennent en M₁ et M₁’. On peut donc décomposer le déplacement de M’en une déformation pure qui amène M’en M₀’ suivie d’un déplacement correspondant à la translation et à la rotation d’ensemble rapportées au point M. Les déplacements élémentaires des particules du corps ayant subit la déformation sont en général exprimées à l’aide de leurs composantes : u, v et w respectivement parallèles aux trois axes de coordonnées x, y et z.

Coordonnées u, v et w des déplacements élémentaires d’un corps

Notion de tenseur des déformations

Le tenseur des déformations peut être traité de la même manière que le tenseur des contraintes. Les règles régissant les contraintes peuvent être transposés aux déformations. Ce tenseur est une matrice carrée 3*3, symétrique noté e_ij qui dans un système d’axes quelconque est représenté par :

Il existe cependant trois directions principales le long desquelles la distorsion est nulle, la déformation est donc entièrement définie par les variations unitaires de longueur dans chacune de ces trois directions.

Comme le tenseur des contraintes, il présente aussi trois invariants J₁, J₂ et J₃ dont les expressions en système d’axes principal sont :

ü J₁ = e₁ + e₂ + e₃

ü J₂ = e₁e₂ + e₂e₃ + e₃e₁

ü J₃ = e₁ e₂ e₃

L’étude des déformations et des contraintes est un aspect essentiel dans la modélisation des milieux continus, principalement pour ses applications à la mécanique des structures déformables notamment dans la formulation des lois de comportement

Lois de comportement mécanique des matériaux

La loi de Hooke

La loi de Hooke va relier les deux tenseurs précédemment définis. Pour étudier ces relations, considérons un parallélépipède rectangle soumis à l’action d’une force normale s_x uniformément répartie sur deux faces opposées.

Action de la force normale sur les faces opposées du parallélépipède

L’expérience prouve que dans le cas où le corps est isotrope ces contraintes normales ne déterminent aucune distorsion des angles du solide élémentaire. Notons d’une part que l’allongement e_x dans le sens des x est lié à la contrainte s_x par la relation :

Dans cette relation E est le module d’élasticité longitudinale ou module d’Young (c’est une constante) qui, pour un matériau homogène isotrope, lie la contrainte à la déformation. Il s’exprime en MPa dans le système international. En effet, il caractérise la raideur de la matière ; à contrainte égale, un matériau ayant un module d’élasticité élevé subira une déformation plus faible qu’un matériau ayant un module d’élasticité petit. Pour déterminer le module E d’un matériau isotrope, on réalise un essai en traction et on enregistre la courbe contrainte = f (déformation). La pente de la courbe dans sa partie linéaire correspond au module élasticité du matériau. D’autre part l’allongement de l’élément dans le sens des x s’accompagne d’une contraction latérale :

e_y = e_z = -n e_x

=

où n est le coefficient de Poisson.

Si on considère par ailleurs que le matériau est soumis à l’action de contraintes normales s_x, s_y et s_z uniformément réparties sur les facettes principales, l’expérience prouve qu’on peut appliquer la méthode de superposition.

Ceci se traduit par les équations donnant l’expression générale de la loi de Hooke pour un corps isotrope. Ces équations reliant la contrainte à la déformation sont :

e_x = 1/E [s_x - n (s_y + s_z)]

e_y = 1/E [s_y - n (s_z + s_x)]

e_z = 1/E [s_z - n (s_x + s_y)]

Ces différentes expressions montrent que la déformation d’un milieu élastique homogène et isotrope ne peut être déterminée, sur la base des contraintes, que si deux modules de déformations sont connus à savoir le module d’Young et le coefficient de poisson.

Comportement élastique, élastoplastique et plastique

Les lois de comportement calées parfaitement sur des essais conduisent à des calculs très lourds et complexes. On utilise ainsi le plus souvent des lois simplifiées reliant contraintes et déformations et qui se prêtent mieux à l’établissement de codes de calcul fournissant ainsi une bonne image du comportement réel du sol.

Les modèles les plus couramment utilisés sont :

Le modèle élastique : c’est le modèle de comportement de matériau le plus simple. Il se caractérise par une relation linéaire (réversible et indépendante de l’histoire des sollicitations) entre la contrainte totale et la déformation (lois de Hooke). Le matériau retrouvera sa forme initiale après relâchement des contraintes externes, à condition toutefois que celles-ci soient inférieures à la limite élastique.

Modèle plastique : à partir d’un premier seuil appelé limite élastique, la déformation augmente à contrainte relativement constante. On note une déformation permanente du matériau.

Modèle plastique parfait

Modèle élastoplastique : en complément de la déformation élastique définie par la loi de Hooke, les modèles élastoplastiques intègrent un degré de déformation plastique permanent. Le matériau présente donc une élasticité linéaire jusqu’à un certain seuil (limite élastique) puis devient parfaitement plastique. La déformation totale se décompose alors en deux parties une élastique et un autre plastique.

Modèle de comportement élastoplastique

On remarquera que les matériaux parfaitement plastiques et élastoplastiques continuent de se déformer même sans accroissement de contraintes.

Les types de comportement d’un matériau

En résumé cette figure donne une illustration de ce qui se passe au delà du comportement élastique : entre les états A et B, le matériau est élastique; il deviendra ensuite plastique si l’on augmente la déformation jusqu’à C. Relâcher la contrainte en C fera revenir le matériau dans une nouvelle position d’équilibre D, selon une pente égale à celle de son comportement élastique.

Critères de plasticité usuels

La grande diversité des matériaux réels se traduit par l’existence d’une multitude de lois de comportement et en particulier d’une grande variété de critères et de lois d’évolution aussi bien en élastoplasticité qu’en élasto-viscoplasticité. Il est illusoire de vouloir établir une liste exhaustive des modèles, d’autant plus que les chercheurs continuent encore à proposer de nouvelles versions. Ce paragraphe sera consacré à une tâche plus modeste qui consiste à présenter le cadre général d’écriture des critères de rupture de Tresca, de Von Mises, de Mohr-Coulomb et de Drucker-Prager.

Critère de rupture de Tresca

Il est d’usage, pour définir mathématiquement le domaine d’élasticité d’introduire une fonction scalaire du tenseur des contraintes, telle que, dans l’espace des contraintes :

F < 0 correspondant à l’intérieur du domaine

F = 0 à sa frontière

F > 0 à l’extérieur du domaine

La fonction F est appelée fonction de charge du matériau dans l’état initial. On désigne aussi couramment par critère de limite d’élasticité, ou critère de plasticité la condition F = 0. Pour le critère de Tresca, la fonction de charge correspondante s’écrit sous la forme :

F = supís_i - s_j - s_0; i, j = 1, 2, 3ý

(s_i, i = 1, 2,3 : contraintes principales)

s₀ apparaît comme la limite d’élasticité en traction simple; elle est égale à l’opposé de la limite en compression simple. Ce critère est aussi appelé « critère en scission maximale », terminologie dont l’origine est évidente si l’on se reporte par exemple à la représentation des contraintes dans le diagramme de Mohr ; sur toute facette l’intensité de la contrainte de cisaillement est bornée par :

où

s₀ / 2 est la résistance en scission simple.

Courbe intrinsèque pour le critère de Tresca

La courbe intrinsèque est par définition la courbe enveloppe de l’ensemble des cercles de Mohr. Dans l’espace des contraintes principales, le domaine d’élasticité du matériau est représenté par un prisme hexagonal régulier.

Domaine d’élasticité du critère de Tresca dans l’espace des contraintes principales

Critère de rupture de Von Mises

Comme le critère de Tresca, il s’agit d’un critère valable pour les matériaux isotropes. Il est également indépendant de la composante sphérique du tenseur des contraintes. La fonction de charge ne dépend donc que du déviateur du tenseur des contraintes. De façon générale pour un matériau isotrope, la fonction de charge qui est fonction des invariants des contraintes s’exprime de manière équivalente par :

Si F est indépendante de la partie sphérique du tenseur des contraintes et ne dépend que de sa partie déviatorique, elle s’exprime en fonction des seules invariantes J₂ et J₃.

Le critère de Von Mises correspond pour un critère de ce type à la forme la plus simple. La fonction de charge s’écrit :

où k est la limite d’élasticité en scission simple et est constant pour un matériau isotrope.

Dans l’espace des contraintes principales (s₁, s₂, s₃), le domaine d’élasticité du matériau est un cylindre droit de rayon .

Domaine d’élasticité dans l’espace des contraintes principales

Il est utile, pour les applications pratiques, de donner l’expression de la fonction de charge de Von Mises, en fonction des contraintes principales :

Cette formule permet d’éviter d’effectuer des calculs fastidieux.

Le critère de rupture de Mohr-Coulomb

Il est apparenté au critère de Tresca, il possède deux propriétés importantes :

1. Il ne porte que sur les contraintes principales extrêmes,

2. Seule la différence entre ces contraintes principales intervient dans l’expression de la fonction de charge.

Ainsi, en ordonnant les contraintes principales s₁, s₂, s₃ sous la forme s₁ ³ s₂ ³ s₃ autrement dit en traction positive, la fonction de charge du critère de Tresca s’écrit sous la forme :

F = s₁ - s₃ - s₀

Le critère de Coulomb conserve la première propriété du critère de Tresca et les deux contraintes principales extrêmes y interviennent explicitement. La fonction de charge de ce critère s’écrit alors sous la forme :

F = s₁ - s₃ + (s₁ + s₃) x sinf - 2ccosf

Elle fait donc intervenir deux constantes caractéristiques du matériau : c (la cohésion) et f (l’angle de frottement interne). Le matériau est dit pulvérulent ou purement frottant si f ¹ 0 et c=0 et il est purement cohérent si f =0 et c ¹ 0 et dans ce cas, le critère se réduit à celui de Tresca avec s₀= 2c.

Le domaine d’élasticité pour ce critère est représenté dans l’espace des contraintes principales par la figure ci-dessous.

Domaine d’élasticité de Mohr-Coulomb dans l’espace des contraintes principale.

Ce critère est sous-tendu par la notion de frottement et suppose que le cisaillement maximal que peut subir le matériau (t) est d’autant plus grand que la contrainte normale en compression est élevée. La limite admissible constitue la courbe intrinsèque dans le plan de Mohr.

Représentation du critère de Mohr-Coulomb dans le plan de Mohr.

Critère de rupture de Drucker-Prager

Au lieu du critère de Coulomb, Drucker et Prager (1952) ont proposé d’utiliser un critère mathématiquement plus régulier qui lui est apparenté et qui s’exprime en fonction simple des invariants I₁ et J₂ du tenseur des contraintes. La fonction de charge pour ce critère est de la forme :

où H est la limite théorique d’élasticité en traction triple (H ³ 0) et a est une constante positive qui dépend du matériau :

Le domaine d’élasticité correspondant est représenté par un cône de révolution au lieu de la pyramide hexagonale du critère de Coulomb

Domaine d’élasticité du critère de Drucker-Prager.

Le modèle de rupture de Cam Clay (modèle du Cambridge)

Le modèle de Cam-clay a été développé par une équipe de chercheurs de l’université de Cambridge. Après avoir posé en 1958 les concepts fondamentaux d’état limite et d’état critique, ces chercheurs ont longuement travaillé à développer un modèle de comportement sur des argiles remaniées isotropes et afin de rendre le modèle quantitativement utilisable, ils ont introduit deux hypothèses importantes qui ont permis de définir la surface d’état limite ainsi que les déformations développées le long d’un cheminement d’état quelconque :

ü La première hypothèse voulant que le sol soit un matériau plastique développant un écrouissage isotrope,

ü La deuxième concernant l’équation de l’énergie dissipée.

Historiquement le nom ‘‘Cam clay’’ a été attribué par Roscoe, Schofied et Thurairajah (1963) à un modèle de comportement, qui sur la base de notion théoriques telles que la surface d’état limite (SEL) et de ‘‘l’indice des vides critiques’’, généralise le comportement de l’argile saturé et humide observé dans les différents essais de laboratoire. Ces auteurs ont analysé l’évolution de la déformation en terme de transfert d’énergie. Ils supposent d’abord que la loi d’écoulement est applicable à l’argile humide, que l’argile humide se comporte comme un milieu isotrope, élastique et non linéaire à l’intérieur de la surface d’état limite (SEL) et qu’il y a des déformations volumétriques réversibles et irréversibles tandis que les déformations de cisaillement sont irréversibles. Ils émettent ensuite une autre hypothèse à savoir la dissipation de l’énergie se fait par frottement interne avec un certain coefficient de frottement

qui est constant ; p et q sont donnés respectivement par les relations :

la contrainte moyenne p peut être donner en terme de contraintes effectives et on notera p’.

est le deuxième invariant déviatorique et peut aussi s’écrire sous la forme :

A la différence des autres modèles, ils introduisent donc un autre espace de représentation qui est l’espace (p’, q, e), e représentant l’indice des vides.

La fonction de charge pour ce critère est donnée dans le plan (p’, q) par la relation :

est la contrainte de préconsolidation.

Burland (1965) a apporté au modèle de ‘‘Cam-clay’’ une modification importante tout en conservant les hypothèses fondamentales de la méthode. On parle ainsi du modèle de ’’Cam-clay modifié’’ dont la fonction de charge s’écrit sur la forme :

La représentation graphique de ces fonctions dans le plan (p’, q) donne la figure sur laquelle on peut remarquer une ligne d’état critique (LEC) correspondant à la droite d’équation q = Mp’. Le diagramme œdométrique de ces fonctions a donné une ligne de consolidation normale (C_c), une ligne de surconsolidations (C_s) les lignes d’état critique correspondant aux lois de plastification.

Cam-clay

Cam-clay modifiée

Les lois de plastification du modèle de Cam-clay dans l’espace p-q’

Diagramme oedométrique des modèles de Cam-clay

La réponse mécanique d’un massif de sol dépend essentiellement des lois de comportement qui lui ont été affectées. Il est donc très important de choisir le modèle rhéologique qui pourra reproduire avec une fidélité optimale la réponse du massif face aux sollicitations qui lui sont imposées et en fonction des grandeurs recherchées. Les modèles de Cambridge sont fondés sur une théorie simple et n’utilisent que quelques paramètres courant de mécanique des sols faciles à déterminer. Le comportement d’un massif de sol, d’après ces modèles, peut être assimilé avec assez de succès à celui d’un matériau élastoplastique dont les relations fondamentales contraintes déformations ainsi que l’équation de la surface d’état limite sont bien définies.

Modélisation par la méthode des éléments finis

Historique et définition de la méthode

Dans son sens le plus général, le terme « modélisation » désigne la conception et l’établissement d’un modèle théorique. Elle consiste à formaliser mathématiquement un processus afin d’être en mesure d’utiliser les avantages de la programmation sur ordinateur. La modélisation s’attache donc à établir un modèle qui est une représentation mathématique simplifiée destinée à mieux appréhender un phénomène complexe, son évolution et les relations et interactions qui le régissent.

La simulation numérique touche actuellement un large éventail de disciplines scientifiques et diverses méthodes ont été développées parmi lesquelles nous pouvons citer :

La méthode des « éléments limites » dans laquelle seule les limites de la zone étudiée sont divisées en éléments, l’intérieur de la masse rocheuse étant représenté mathématiquement comme un milieu continu infini ou fini suivant que la frontière est ouverte ou fermée. Cette méthode est souvent utilisée en mécanique des roches.

La méthode «des éléments finis» (MEF ou FEM) qui est la méthode la plus efficace et la plus utilisée. Elle consiste à utiliser une approximation simple de la géométrie et des variables décrivant le phénomène physique telles que les déplacements, la vitesse, la pression, la température, afin de ramener le problème continu comportant une infinité d’inconnues vers un système algébrique facile à résoudre et à nombre fini d’inconnues ;

La méthode des «éléments discrets»

C’est une méthode générale qui s’applique aux problèmes de mécanique des milieux continus comprenant :

q Les problèmes de mécanique des solides (élasticité, plasticité, mise en forme…)

q Les problèmes de mécanique des fluides, d’hydraulique.

q Les problèmes d’électromagnétisme.

L’origine de cette méthode est à rechercher dans les études des structures en « treillis » communément faites par les techniciens du génie civil. Ces structures en treillis sont composées de parties indépendantes reliées entre elles à des nœuds. L’analyse de telles structures s’effectue en considérant le comportement de chaque partie indépendante, puis en assemblant ces différentes parties de telle sorte que l’équilibre des forces et la compatibilité des déplacements soit satisfaite en chaque nœud. D’une manière beaucoup plus générale, nous pouvons dire que la MEF est un outil de discrétisation. Le domaine dans lequel on désire effectuer la simulation est découpé en « élément » dont l’ensemble constitue le maillage de la structure. On dit alors que la structure est maillée ou discrétisée. L’application des conditions d’équilibres et des lois de comportement de la mécanique permet de construire des équations approchées dont les inconnus sont les valeurs de la solution en un ensemble bien choisi de points appelés nœuds de la discrétisation. Ainsi à la différence de la mécanique des milieux continus (théories de l’élasticité), qui utilise des fonctions représentant les divers grandeurs (contraintes, déformations, déplacements) en tous points de la structure, le calcul par éléments finis s’applique à un nombre fini de grandeurs que l’on traite matriciellement. Dés lors, il est indispensable de substituer à la structure étudiée, un modèle mathématique composé d’un certain nombre éléments, appelés éléments finis, reliés entre eux en un nombre fini de point (nœuds). Les éléments obtenus par maillage restent des corps déformables continus mais leur forme plus simple permet de les étudier plus aisément.

Lorsque la discrétisation est terminée, le stade suivant consiste à calculer les caractéristiques des éléments, c’est à dire les relations entre les forces et les déplacements généralisés aux nœuds par l’intermédiaire de la matrice de rigidité qui établit une relation entre le vecteur des forces et celui des déplacements au nœud considéré. Cette relation est de la forme :

où est le vecteur des efforts généralisées, celui des déplacements correspondants et la matrice de rigidité.

La matrice de rigidité est une matrice singulière cela signifie mathématiquement qu’elle est symétrique, que les termes diagonaux sont positifs et que la somme des termes de chaque colonne est nulle.

Une simulation réaliste peut exiger des centaines de milliers de nœuds et d’éléments. Les systèmes à résoudre sont donc énormes et il va de soi qu’un tel traitement requiert des calculs considérables dont seul l’outil informatique permet de venir à bout. Ce type de programme comporte grossièrement les parties suivantes :

ü Introduction des données :

· Géométrie de la structure

· Les matériaux qui la composent

· Conditions de chargement

· Conditions aux limites

ü Bibliothèque des éléments : elle comprend un certain nombre d’éléments et leur modèle mathématique (leur matrice de rigidité). Ce sous programme permet de calculer la matrice de rigidité en fonction des données (coordonnées nodales et caractéristiques des matériaux).

ü La résolution du système : permet d’obtenir les contraintes, les déformations et les déplacements nodaux.

5.2. - La démarche de la méthode des éléments finis

La démarche de la méthode peut être exposée autrement. Pour ce faire, on considère un domaineW. La méthode consiste à rechercher une solution approchée de la solution exacte sous forme de champ défini par morceaux sur des sous domaine deW. Les n domaines W_i doivent être tels que :

Et

" i ¹ j

où désigne l’intérieur de W_i. Autrement dit, les W_i sont une partition deW.

Le champ dans chaque sous domaine W_i est déterminé par un nombre fini de valeur du champ (ou de valeur de ses dérivées) en des points choisis arbitrairement dans le sous domaine, et appelés nœuds. Le champ local est une interpolation entre les valeurs aux nœuds. Chercher une solution par éléments finis consiste donc à déterminer quel champ local on attribue à chaque sous domaine pour que le champ global obtenu par juxtaposition de ces champs locaux soit proche de la solution du problème.

Solution approchée discontinue

bolutions approchées continues

Méthode d’approximation d’une solution réelle

La qualité de la solution approchée dépend de la division en sous domaines (nombre et dimensions des sous domaines), du choix de la famille de champs locaux dans chaque sous domaines, et des conditions de continuité qu’on impose aux frontières des sous domaines. Une fois ces choix faits, il reste à rechercher, une combinaison de champs locaux qui satisfait approximativement les équations. Pour résoudre un problème par la méthode des éléments finis, on procède par étapes successives qui peuvent être schématisées par l’organigramme suivant :

Organigramme d’un programme d’éléments finis

Aujourd’hui, le travail à effectuer dans ces différentes étapes est assisté par les logiciels. Ainsi de très nombreux programmes d’analyse de structure par éléments finis ont été mis au point ces dernières années et à côté d’un nombre impressionnant de petits programmes simples destinés à des applications spécifiques limitées, il existe quelques systèmes complexes destinés à un emploi généralisé. Parmi ces derniers citons : PLAXIX 7.2^© MAC STAR 2.2^©; NASTRAN (NASA, USA), CASTEM (CEA), etc. Ils se présentent tous avec le même organigramme. Il reste que pour maîtriser leur utilisation, il est indispensable de comprendre les fondements de la méthode afin de choisir intelligemment parmi les options qu’ils proposent. Le choix des logiciels est guidé par un certain nombre de critères parmi lesquels on a :

• La documentation : note théorique décrivant les possibilités du programme, les modèles choisis, les formulations et algorithmes et en fin la disponibilité de la notice d’utilisation (qui doit être bien rédigé pour évité d’éventuel perte de temps) ;
• Les possibilités du logiciel : type d’analyse possible, sauvegarde et restauration des données ;
• Autres critères : structure déjà traitées, entretien, aide à l’utilisation, etc.

Pour clore ce chapitre nous disons qu’une des forces de la méthode des éléments finis est que le choix initial des éléments peut être revu à volonté en particulier pour raffiner le calcul sur une partie du système sur laquelle on choisira un maillage plus serré. Cependant la validité des résultats que cette méthode fournit dépend de la pertinence des choix réalisés à chaque étape de la modélisation et de la mise en œuvre. Il est important aussi d’ajouter que la simulation numérique, sans totalement se substituer aux essais expérimentaux, constitue un outil complémentaire qui s’avère souvent beaucoup plus rapide et moins onéreux que ces derniers.

Par la méthode des éléments finis, on peut obtenir des calculs précis des contraintes et des déformations dans un massif de sol soumis à des charges uniformes ou complexes à condition d’en connaître la loi de comportement et les conditions aux limites. Cette méthode peut donc être appliquée à des matériaux ayant un comportement élastoplastique, comportement dont les équations de la surface d’état limite ainsi que les relations contraintes déformations sont clairement définies.