PARTIE 1 - ÉLÉMENTS THÉORIQUES

$translationBooks

Chapitre A. Généralités

Chapitre A. Généralités

Lors du calcul d’une structure par éléments finis, le cadre des hypothèses de modélisation doit être posé ainsi que les objectifs à atteindre, qui sont généralement la détermination de :

Parmi les hypothèses de modélisation, il est préférable de traiter, par ordre chronologique :

A1. Formulation générale pour les calculs élastiques linéaires

A1. Formulation générale pour les calculs élastiques linéaires

A2. Dimensionnalité de la modélisation

A2. Dimensionnalité de la modélisation

A3. Choix des éléments finis

A3. Choix des éléments finis

A4. Interaction entre la structure et son environnement

A4. Interaction entre la structure et son environnement

A5. Estimation de la qualité de la solution numérique approchée

A5. Estimation de la qualité de la solution numérique approchée

$translationBooks

A1. Formulation générale pour les calculs élastiques linéaires

A1. Formulation générale pour les calculs élastiques linéaires

La formulation qui n’est plus réellement à faire car la majorité des logiciels proposés de nos jours s’appuie sur une formulation en déplacement car la prise en compte des conditions aux limites ne pose pas de difficulté notoire à l’inverse d’une formulation en contrainte. Il existe également la formulation mixte en déplacement et en contrainte. Ces trois formulations sont les principales, mais il existe d’autres formulations que nous n’évoquerons pas dans ce chapitre.

Dans ce document, nous ne présentons que la formulation en déplacement (la plus utilisée) : dans ce cadre, les inconnues de base sont les déplacements aux nœuds, notés classiquement q, dont sont déduits tous les autres résultats (déformations ε, contraintes σ …).

Les déplacements en tout point d'un élément, notés u(x,y,z), sont obtenus par interpolation des valeurs aux nœuds, à l'aide des fonctions d'interpolation (appelées fonction de forme), notées N(x,y,z), polynomiales de degré faible (1 à 3 en général) :

La relation déformations - déplacements s'obtient par dérivation à l'aide de l'opérateur (matrice) de dérivation D (différent selon le contexte), faisant apparaître la matrice B reliant déformations et déplacements aux nœuds :

Pour obtenir les contraintes, il est enfin indispensable d'introduire la loi de comportement du matériau constitutif. Ainsi, en élasticité linéaire, la relation contrainte – déformation s'écrit (avec H la matrice de Hooke) :

On constate donc que toutes les quantités d'intérêt (déplacements en tout point, déformations, contraintes) s'obtiennent à partir de la connaissance des déplacements aux nœuds du maillage

Remarque : comme les déformations (et donc les contraintes) sont obtenues par dérivation des déplacements aux nœuds, il y a dégradation de la précision lors du calcul de ces grandeurs.

Pour un calcul statique en mécanique des structures, la détermination de la solution numérique du problème par la Méthode des Éléments Finis (MEF) peut se résumer en 4 étapes principales :

Étape 1. Détermination des matrices et des vecteurs élémentaires sur chaque élément de volume V par les relations suivantes (la notation BT désignant la transposée de la matrice B) :

 matrice élémentaire de rigidité

 vecteur élémentaire des charges équivalentes

avec  le vecteur des charges volumiques et  celui des charges surfaciques. Le vecteur des charges équivalentes permet de « ramener aux nœuds » les charges appliquées aux éléments (dans le volume ou en surface), selon l'interpolation retenue

Attention, cela ne correspond en général pas à l'équirépartition aux nœuds de la charge totale appliquée à un élément.

Étape 2. Calcul de la matrice de rigidité K et du vecteur des charges F à l’échelle de la structure par assemblage des quantités élémentaires et introduction des conditions aux limites.

Étape 3. Détermination du vecteur q des déplacements nodaux par la résolution du système linéaire (de grande taille) :

Le vecteur q donne la valeur du déplacement aux nœuds du maillage, comme c'est le cas du vecteur des charges équivalentes F.

Etape 4. Déduction des quantités d’intérêt telles que déformations et contraintes par un post-traitement du vecteur q des déplacements nodaux.

En revenant au niveau des éléments, les valeurs des déplacements en tout point, des déformations et des contraintes sont déduites des déplacements aux nœuds à l'aide des relations précédentes. Contrairement aux déplacements, les contraintes et les déformations sont des grandeurs calculées élément par élément (et il n'y a pas toujours continuité de ces quantités d'un élément à l'autre).

Remarque : il est possible de faire une analogie entre résolution par éléments finis d'un problème mécanique et d'un problème thermique en régime permanent linéaire : dans ce cas, les températures aux nœuds T sont obtenues en résolvant le système linéaire : Λ.T = Φ, Λ étant la matrice de conductivité et Φ le vecteur des flux aux nœuds.

Le degré des fonctions de forme N(x,y,z), qui se déduit pour les éléments massifs (ou de mécanique des milieux continus - MMC) du nombre de nœuds sur l'arête d'un élément fini (2 nœuds : éléments linéaires, 3 nœuds : éléments quadratiques) se répercute sur le degré des déformations et des contraintes évaluées par éléments finis. Par exemple, si l'opérateur de dérivation D correspond à des dérivées premières (cas de la MMC), les déformations et contraintes seront constantes par éléments pour des éléments linéaires, linéaires pour des éléments quadratiques… Pour mémoire, en RdM, les déformations font intervenir des dérivées d'ordre 2 (flexion des poutres, plaques…).

image029.jpg

Fonctions de forme - Triangle à 3 noeuds (fonctions linéaires)

image030.jpg

Fonctions de forme - Triangle à 6 noeuds (fonctions quadratiques)

Attention : il n'est pas possible de mélanger dans un même maillage des éléments finis dont les fonctions de forme sont de degrés différents (par exemple éléments linéaires et quadratiques). En effet, la continuité des déplacements ne peut être assurée à l'interface d'éléments de fonctions de forme de degrés différents.

Dans les logiciels, la plupart des éléments finis sont dits isoparamétriques, c'est-à-dire que les différents éléments finis d'un maillage sont construits par transformation géométrique d'un élément de référence (lui-même exprimé dans un repère de coordonnées unitaires ξ,η,χ).

image031.jpg

Élément isoparamétrique (quadratique) : différentes transformations (2D)

Ces transformations géométriques, reposant en général sur les mêmes nœuds que l'interpolation sur les déplacements (d'où le qualificatif « isoparamétrique » : mêmes paramètres), sont en particulier caractérisées par leur matrice jacobienne J (matrice des dérivées partielles du premier ordre de la transformation géométrique).

image032.jpg

Élément isoparamétrique (quadratique) : repère local et global (2D)

Un des avantages de cette transformation est la possibilité d'approcher les frontières du domaine d'étude par des géométries polynomiales. Cependant, l'inconvénient principal est de complexifier le calcul des matrices et vecteurs élémentaires. Leur évaluation se fait alors par intégration numérique à l'aide d'un schéma de type quadrature : somme pondérée de valeurs en certains points situés à l'intérieur des éléments (poids ω_i, points ξ_i) caractéristiques de la méthode (méthode de Gauss, de Hammer…). Dans ce cas, le calcul de la matrice de rigidité élémentaire et du vecteur élémentaire des charges équivalentes (aux charges volumiques) s'écrivent :

   

Quand les matrices (faisant intervenir la matrice B) sont évaluées de cette façon, il est naturel de calculer aux mêmes points les contraintes et les déformations (également fonction de B). Les logiciels fournissent alors ces grandeurs aux points d’intégration (aux points de Gauss ou de Hammer selon le schéma d'intégration) ; de plus, parler de valeurs aux nœuds d'un maillage pour les contraintes ou les déformations n'a pas de sens, car il y a autant de valeurs en un nœud qu'il y a d'éléments concourants en ce point et il n'y a pas forcément continuité d'un élément à l'autre ! Vouloir à tout prix obtenir une valeur (de contrainte ou de déformation) en un nœud du maillage suppose donc un traitement : extrapolation à partir des valeurs aux points d'intégration et moyenne des valeurs provenant des différents éléments concourants, qui conduit au minimum à un lissage des résultats obtenus (cependant, parler de valeur aux nœuds d'un élément particulier est justifié). Quand ces quantités sont obtenues par intégration analytique (cas des poutres par exemples), elles peuvent être exprimées aux nœuds, au centre de gravité… les remarques précédentes sur la continuité entre éléments restant valables.

Remarque : certains éléments finis présents dans les codes sont dits sous-intégrés (ou à intégration réduite) : il s'agit d'éléments pour lesquels l'intégration numérique lors du calcul des quantités élémentaires (matrice de rigidité et vecteur des charges équivalentes) est effectuée avec moins de points que celui mathématiquement nécessaire ; cela peut améliorer dans certains cas la précision des calculs (en évitant les phénomènes de bloquage mais également conduire à des problèmes numériques (matrice singulière par exemple ou apparition de modes d’hourglass).

image033.jpg

Éléments quadratiques et points d'intégration (points de Gauss)

La façon de prendre en compte les conditions aux limites en déplacement (encastrement, appuis…) peut entraîner des différences de résultats selon la technique numérique retenue, qui peuvent être résumées à deux : la méthode de pénalisation, simple à mettre en œuvre mais qui présente l'inconvénient d'être sensible à l'ordre de grandeur des termes de la matrice de rigidité et la méthode des multiplicateurs de Lagrange, insensible aux problèmes de précision évoqués pour la première mais augmentant légèrement la taille du système à résoudre.

Remarque : en dynamique, le déplacement, en tout point de coordonnées x, y, z et à tout instant t s'obtient, comme en statique, par interpolation des valeurs aux nœuds : u(x,y,z,t)=N(x,y,z).q(t), elles-mêmes vérifiant l'équilibre qui s'écrit, à chaque instant t :

la matrice de masse étant obtenue par assemblage des matrices de masse élémentaires  (masse volumique) et C étant la matrice d'amortissement (visqueux). Ceci nécessite la résolution d'un système d'équations différentielles par des méthodes numériques adéquates (cf. chapitre dynamique).

$translationBooks

A2. Dimensionalité de la modélisation

A2. Dimensionnalité de la modélisation

Il est indispensable de transformer le problème réel autrement dit de le simplifier pour modéliser l’interaction entre la structure et son environnement. Pour ce faire, il existe 2 techniques :

  • La première est de passer d’un problème spatial à un problème de dimension réduite tel que l’axisymétrie ou le 2D plan (en contraintes planes ou en déformations planes). Elle permet une réduction de l’espace (cf. §2.2).
  • La seconde permet de rester dans la dimension 3D mais en tenant compte d’une réduction du modèle grâce à des hypothèses cinématiques. Ceci est constaté lors du passage à la théorie des poutres, des plaques ou des coques (cf. §2.1).

Ces 2 techniques induisent des coûts de calcul réduits. En revanche, l’utilisateur doit faire très attention à la simplification employée car cette dernière fait appel à des hypothèses qui doivent rester dans le domaine de validité du problème réel afin d’obtenir des résultats pertinents.

1) Cas des éléments finis de RDM

Du point de vue des éléments finis, la différence essentielle entre MMC et RdM concerne la géométrie qui est simplifiée moyennant des hypothèses supplémentaires : le problème tridimensionnel initial est alors ramené à un problème bidimensionnel (surface moyenne pour les plaques et coques) ou unidimensionnel (fibre moyenne pour les barres et poutres), mais représenté dans un espace tridimensionnel (à la différence des problèmes plan, cf. ci-dessous). Lorsque la MEF est utilisée pour résoudre un problème de RdM, les éléments finis sont donc des éléments spécifiques, pour lesquels il faudra fournir des caractéristiques géométriques (section, inerties pour les poutres, épaisseur pour les plaques et coques). De plus, ils combinent les effets de traction/compression (pour les poutres), ou de membrane (pour les plaques et coques), qui sont traités de manière similaire à la MMC, aux effets de flexion, qui sont traités spécifiquement.

Pour la partie flexion, la simplification géométrique induit une définition particulière des déformations, ce qui entraîne pratiquement des expressions différentes de l'opérateur de dérivation D selon les cas. Une autre conséquence porte sur la définition des inconnues aux nœuds : en effet, alors qu'en MMC, les inconnues sont les composantes du déplacement, en RdM, on est amené à rajouter des inconnues correspondant aux rotations, puisqu'il n'est plus possible de les évaluer directement à partir des déplacements aux nœuds en raison de la simplification géométrique. Le choix initial de traiter le problème posé en MMC ou RdM implique donc le choix d'éléments finis dans des familles différentes ; il est de plus a priori interdit de les mélanger car, à l'interface entre éléments de natures différentes, on se retrouve avec des rotations non transmises, sauf disposition spécifique. De plus, les contraintes calculées sur ces éléments sont en général des contraintes généralisées (ou efforts de la RdM : effort normal, tranchant, moment de torsion, fléchissant). Pour obtenir les contraintes (de la MMC) en un point, il est nécessaire de fournir des informations supplémentaires (position dans la section de la poutre par exemple).

Dans ce cadre, le choix d'un élément fini présente une difficulté supplémentaire liée à la prise en compte ou pas de l'énergie de cisaillement (éléments de poutre de Euler-Bernoulli ou de Timoschenko, éléments de plaque de Love-Kirchoff ou de Reissner-Mindlin), ce choix étant lié à des considérations géométriques (élancement de la section de la poutre ou épaisseur de la plaque). De plus, dans le cas où on choisit de prendre en compte l'énergie de cisaillement, des problèmes numériques peuvent intervenir (blocage en cisaillement), qui rendent certains éléments finis d'utilisation délicate


L'élément fini de poutre d'Euler-Bernoulli permet de représenter exactement un moment fléchissant variant linéairement le long de la fibre moyenne d'un élément (les fonctions de forme étant de degré 3 pour les déplacements de flexion et le moment obtenu par dérivée seconde des déplacements) : il n'est donc pas utile de mettre plusieurs éléments entre nœuds chargés par des forces ponctuelles (alors qu'il faut mailler finement lorsque le chargement est réparti entre 2 points).


Enfin, pour les éléments de plaque et coque, la convergence monotone n'est pas toujours assurée selon la forme du maillage, ce comportement étant lié à la formulation même des éléments. Ce type de comportement est illustré sur la figure ci-dessous, montrant l'évolution de l'erreur relative sur la flèche ω et le moment Mx en fonction du logarithme du nombre de degrés de liberté (dans le cas d'une plaque rectangulaire appuyée sur ses quatre côtés soumise à une charge uniforme) : en vert, un élément non-conforme (COQ3) dont les résultats sont bien moins précis qu'un élément conforme (DKT), et de plus dont la convergence n'est pas monotone (pour les moments Mx sur la figure) : plus d'éléments conduit paradoxalement à un résultat qui peut être moins précis !

Plaque en flexion : convergence (flèche ω et moment Mx) en fonction du log du nombre de degrés de liberté

Ce domaine, particulièrement important en génie civil, présente donc des difficultés spécifiques qui seront développées dans les prochains chapitres.

2) Calculs bidimensionnels

Les problèmes étudiés sont par nature tridimensionnels ; cependant, il est plus rapide d'effectuer des calculs bidimensionnels. Dans certains cas, il est possible de ramener l'étude d'un problème tridimensionnel à celui d'un problème bidimensionnel :

  • si le problème admet un axe de révolution (pour la géométrie, le chargement et les conditions aux limites) : il est possible de faire un calcul axisymétrique pour lequel aucune hypothèse supplémentaire n'est faite ; dans le cas où le chargement n'est pas axisymétrique, il est possible de le décomposer en séries de Fourier et de traiter le problème initial par superposition de calculs axisymétrique ;

image050.jpg 
Cylindre sous pression ➡ axisymétrie

  • si on fait l'hypothèse de négliger les contraintes ou les déformations hors plan selon que la structure a une épaisseur très faibles ou très importante : on se place donc en hypothèse de contraintes planes ou de déformations planes respectivement ; la solution obtenue est alors une approximation du problème tridimensionnel. Il faut garder en tête qu'en contraintes planes, les déformations hors plan sont non nulles (idem en déformations planes pour les contraintes hors plan).

image051.jpg Barrage : épaisseur importante ➡ déformations planes

 Assemblage : faible épaisseur ➡ contraintes planes

3) Prise en compte de symétries

Certains problèmes présentent des symétries (axe de symétrie, plan de symétrie…) et il peut être intéressant d'en profiter pour rendre les calculs par éléments finis plus rapides. Il faut cependant se souvenir que la résolution ne fournira que des solutions elles-mêmes symétriques (notamment lors de calculs de modes propres).

Pour cela, il est important d'avoir en tête que, comme mentionné ci-dessus, la symétrie doit concerner aussi bien la géométrie, que le chargement et les conditions aux limites pour pouvoir être prise en compte dans le calcul. La solution de la partie de la structure modélisée ne représentera cependant la solution du problème complet que si les conditions de symétrie sont intégrées au modèle. En mécanique des solides (MMC), celles-ci consistent à imposer que les composantes du déplacement perpendiculaires à l'axe/au plan de symétrie sont nulles. En RdM, il convient de penser à ajouter la condition de nullité des rotations autour de l'axe perpendiculaire à l'axe/au plan de symétrie. Enfin, la prise en compte de la symétrie peut avoir des conséquences sur le chargement : par exemple, il faut penser à n'appliquer que la moitié de l'intensité d'une charge ponctuelle appliquée sur un axe de symétrie.

$translationBooks

A3. Choix des éléments finis

A3. Choix des éléments finis

Le choix de l’élément est une étape importante. L’objectif est de sélectionner le type d’élément (à savoir sa forme et l’ordre des fonctions de forme associées) et la taille de l’élément. Le type et la taille des éléments définissent l’allure et la précision des champs de déplacement, donc, de déformation et de contraintes.

En plus, de la forme de l’élément, vient s’ajouter son aspect. Il faut éviter de produire des éléments dégénérés ou trop dénaturés (aplatis ou allongés) car cela dégrade la précision de la résolution du problème.

De manière plus globale, le maillage généré peut être structuré (répartition régulière des éléments) ou non structuré. Il est possible de combiner sur un même domaine des parties structurées et non structurées suivant la complexité géométrique.

Dans le même esprit, la taille de l’élément choisi est dépendante de la géométrie de la structure à mailler et des chargements qui lui sont appliqués. Les zones à fortes variations (gradients) de contrainte ou à contraintes élevées (contact-frottant ou fissuration par exemple) déterminent les parties nécessitant une finesse de maillage supérieure aux autres parties pour bien percevoir les contraintes ainsi que les déformations.

Une des premières règles à respecter est de lancer une première simulation avec un premier maillage afin de déterminer les zones sensibles énumérées précédemment et ensuite raffiner le maillage où cela s’avère nécessaire.

Le modélisateur doit avoir un regard critique sur son maillage vis-à-vis de la géométrie de la structure et des zones importantes à observer.

Pour bien choisir les éléments finis à utiliser lors du maillage, il est nécessaire de réfléchir au type de calcul souhaité :

  • l'ensemble de la géométrie va-t-elle être représentée, auquel cas le problème relève de la MMC et les éléments sont donc de type « massif » ; dans ce cas, il faut savoir si le problème est tridimensionnel ou bidimensionnel (contraintes planes, déformations planes, axisymétrie… cf. 3.) et si on souhaite utiliser des éléments finis à fonctions de forme linéaires ou quadratiques : la précision des calculs est meilleure avec un maillage d'éléments quadratiques mais a un coût de calcul plus important : il faut donc faire un compromis ; dans tous les cas, les degrés de liberté sont les composantes du déplacement (u, v en 2D, u, v, w en 3D) ; les principaux éléments sont listés dans le tableau ci-dessous ;

  • ou la géométrie est simplifiée, cas de la RdM (ou calcul des structures cf. 2.) : dans ce dernier cas, les éléments finis sont donc
    • des éléments de barre/poutre (une barre ne reprend que de la traction ou de la compression à la différence d'une poutre qui reprend en plus de la flexion ; attention, certains logiciels ne différencient pas les deux et parlent d'éléments de barre pour désigner également des éléments de poutre)
    • des éléments de plaque/coque (la différence entre plaque et coque est liée à la courbure du feuillet moyen et la plupart des logiciels ne font pas la distinction).


Pour ces éléments, en plus des degrés de liberté de déplacement, les éléments de RdM comportent des degrés de liberté de rotation (θx, θy, θz), permettant de prendre en compte la géométrie non maillée (section pour les poutres, épaisseur pour les plaques et coques). De plus, se pose la question de savoir si l'énergie de cisaillement doit être prise en compte ou pas (éléments de poutre de Navier-Bernoulli ou de Timoshenko, éléments de coque de Love-Kirchhoff ou de Mindlin-Reissner). Dans le cas des éléments de coque, comme mentionné auparavant, se pose enfin la question de la qualité des éléments (éléments conformes ou pas). Il est donc particulièrement difficile de faire le choix d'un élément fini de plaque/coque, en particulier lorsque la documentation est succincte ; il peut être judicieux dans ce cas, de réaliser un calcul sur un cas dont la solution est connue afin de tester la qualité des éléments disponibles.

En RdM, les différents éléments généralement rencontrés sont décrits dans le tableau II ci-dessous.

$translationBooks

A.4. Interaction entre la structure et son environnement

A.4. Interaction entre la structure et son environnement

La prise en compte de l’interaction entre la structure et son environnement : cette considération permet de diminuer au mieux les écarts entre la modélisation par éléments finis et la réalité. L’effet de l’environnement sur la structure est déterminé à partir de grandeurs nodales telles que les déplacements nodaux et les efforts nodaux. Les premiers sont liés aux conditions aux limites et les seconds aux chargements extérieurs [1].

  • Les conditions aux limites en déplacement (nodaux) permettent d’imposer aux nœuds une valeur de leur déplacement (nulle ou non nulle). Les déplacements imposés sont souvent appelés contraintes cinématiques. Elles permettent également de lier les déplacements de certains nœuds. Dans un premier temps, il est préférable de déterminer s’il faut prendre en compte la ou les symétries avant d’imposer les conditions aux limites. L’utilisation de symétrie impose de bloquer les déplacements nodaux perpendiculaires au plan de symétrie en 3D ou à l’axe de symétrie en 2D en imposant une valeur nulle (Figure 3) (cf. §4).

Conditions aux limites correspondant à une symétrie par rapport à deux axes dans le cas de la plaque trouée sollicitée en traction

Dans un 2ème temps, il est indispensable de neutraliser les mouvements de corps rigide. Une modélisation par éléments finis correctement définie doit interdire les translations et les rotations libres. En 3D, il faut donc éviter les 6 mouvements de corps rigide énoncés précédemment (en 2D, il s’agit de 3 mouvements, Figure 4).

Mouvements de corps possibles en 2D : a) Translation horizontale possible ; b) Translation verticale possible ; c) Rotation autour de la rotule ; d) Tout mouvement de corps rigide neutralisé

Une fois ces deux étapes réalisées, il est nécessaire de vérifier que tous les mouvements de corps rigide sont correctement annulés et qu’aucun mouvement de corps rigide n’a été bloqué alors que ce mode avait déjà été supprimé. Dans le premier cas, aucun calcul n’est possible et dans le second cas, des déformations non attendues risquent d’apparaître.

  • Le chargement correspond à des efforts extérieurs exercés sur certaines parties du maillage. Parmi les chargements, il existe les forces à distance telles que la pesanteur et les inerties. Elles sont généralement modélisées par des forces volumiques traduites par des forces nodales sur l’ensemble des nœuds du domaine. Ce type de force est souvent simple à mettre en place.

Il existe également les forces de contact telles que les pressions ou toute autre force nécessitant un contact avec la structure. Ces forces peuvent être surfaciques, linéiques ou ponctuelles. L’application de ces forces se traduit également par des forces nodales. Une attention particulière doit être portée à la traduction de ces forces de contact lors de la modélisation vis-à-vis de son domaine de validité. L’utilisation d’une force ponctuelle peut gérer des singularités telles qu’une concentration de contraintes dans le voisinage proche du nœud d’application de la force ponctuelle. Donc pour éviter ce type de singularité, il est nécessaire de répartir la charge par l’intermédiaire d’un patin sur lequel la force ponctuelle sera appliquée. Cela revient à appliquer une force surfacique sur la structure englobant un voisinage plus ou moins grand du point d’application de la force ponctuelle. Se pose alors la question du raffinement de maillage dans cette zone et son influence que les résultats obtenus.

Les conditions de raccords

Il existe différents types d’éléments finis tels que les éléments volumiques, plaques, coques, poutres et barres. Les plaques et coques minces sont des éléments dont l’épaisseur est plus petite que les 2 autres dimensions. Une plaque ne travaille que perpendiculairement à son plan (3 degrés de liberté (ddl) par nœud : une translation et 2 rotations) alors qu’une coque travaille dans son plan et perpendiculairement à celui-ci (6 ddl : 3 translations et 3 rotations). Les éléments volumiques ont 6 ddl par nœud. Les éléments poutres ont 6 ddl par nœud également en 3D. Le chapitre 2 présente des développements pratiques sur ce point.

$translationBooks

A.5. Estimation de la qualité de la solution numérique approchée

A.5. Estimation de la qualité de la solution numérique approchée

L’écart entre la solution exacte du problème et la solution approchée obtenue par la méthode des éléments finis permet de connaître la qualité de la solution : c'est l'« erreur de discrétisation ». La solution exacte n’étant généralement pas connue, l’idée consiste à estimer cet écart en calculant un « estimateur d’erreur ». On distingue :

  • des estimateurs d’erreur globale afin d’évaluer la qualité de la solution sur l’ensemble du domaine (méthode de lissage des contraintes, des résidus d'équilibre, de l'erreur en relation de comportement…),
  • des estimateurs d’erreur locale afin d’évaluer la qualité d’une quantité d’intérêt telle que le déplacement en un point ou la contrainte dans une zone (méthode des résidus d'équilibre pondérés…).

Ces différents outils d’estimateur d’erreur, disponibles selon les logiciels utilisés, pour un coût numérique variable, peuvent être utilisés dans deux buts :

  • améliorer la qualité des résultats d’un calcul par éléments finis en raffinant de façon automatique le maillage et/ou la discrétisation temporelle,
  • obtenir un intervalle de confiance (bornes inférieure et supérieure) associée au calcul de l’erreur globale ou sur une quantité d’intérêt.

Il est important de souligner que certains estimateurs d’erreur, tels que les méthodes de lissage des contraintes, permettent uniquement d’avoir une indication de l’erreur commise alors que d'autres, tels que les méthodes des résidus d’équilibre, permettent de garantir les résultats numériques obtenus via le calcul d’un encadrement de l’erreur.

$translationBooks

Chapitre B. Dynamique

Chapitre B. Dynamique

Pour un certain nombre d’applications, tel que les calculs sismiques, les impacts, les études vibratoires, … il est nécessaire de considérer les phénomènes dynamiques.

Les chargements dynamiques appliqués à une structure de génie civil appartiennent basiquement à deux catégories :

  • les phénomènes assimilables à des phénomènes stationnaires : écoulement permanent de vent, houle, machine tournante,
  • les phénomènes transitoires : impact, explosion, séisme.

Concernant les mouvements sismiques, s’ils sont théoriquement considérés comme transitoire, il est néanmoins admis de les assimiler comme des phénomènes stationnaires pendant leur durée de phase forte. Pour les cas où l'on cherche à modéliser la structure en intégrant des non linéarités géométriques ou matérielles, on ne peut plus considérer de caractère stationnaire

On distingue ensuite les moyens de représenter les catégories de chargements :

  • Stationnaire :
    • Transformée de Fourier complexe (TF) ;
    • Densité spectrale de puissance (DSP) ;
    • Spectre de réponse d’oscillateur (SRO).
  • Transitoire :
    • Courbe de chargement de déplacement, vitesse ou accélération exprimée en fonction du temps ;
    • Effort ou pression exprimé en fonction du temps.

Deux grandes familles d’analyses peuvent être considérées :

  • L’analyse modale qui permet de connaître les fréquences propres et les modes propres d’une structure et servira pour caractériser :
    • la réponse au chargement stationnaire appliqué via une méthode de réponse spectrale ;
    • la réponse temporelle par intégration de Duhamel de chaque réponse modale à la courbe de chargement ;
    • une fonction de transfert convoluée au signal exprimé de manière fréquentielle pour délivrer une réponse en DSP ou TF.
  • La dynamique temporelle qui permet de calculer la réponse dynamique transitoire de la structure pour une excitation temporelle quelconque. Cette résolution est faite à l’aide de schémas d’intégration temporelle, qui peuvent être explicites ou implicites.

Les schémas explicites imposent de choisir des pas de temps très petits ; ils sont donc le plus souvent utilisés pour résoudre des problèmes sur des temps courts (type impact). Au contraire, les schémas implicites permettent d’utiliser des pas de temps plus grands et sont donc privilégiés pour étudier des plages temporelles plus grandes.

Exemples d’applications


Applications Représentation du chargement Grandeurs accessibles
Modal Analyse vibratoire TF TF
DSP DSP
Suivi de fréquence propre SRO Extrema spectraux de quantités d'intérêt variées
Transitoire implicite Etude sismique Accélérations, vitesses, forces, pressions ou déplacements en fonction du temps Quantités d’intérêts diverses exprimées au cours du temps
Ebranlement
Transitoire explicite Chute d’un objet Modélisation de projectiles en contact, chocs Quantités d’intérêts diverses exprimées au cours du temps
Impact d’avion Accélérations, vitesses, forces, pressions ou déplacements en fonction du temps

Le problème dynamique une fois discrétisé par éléments finis se ramène à la résolution de l’équation d’équilibre suivante (cf. chapitre 1) :

  • M la matrice de masse exprimée aux nœuds,
  • C la matrice d’amortissement exprimée aux nœuds,
  • K la matrice de raideur exprimée aux nœuds,
  • q le vecteur des déplacements nodaux,
  • q' le vecteur des vitesses nodales,
  • q¨ le vecteur des accélérations nodales.

Dans le cas de l’analyse modale, on a recours au calcul des pulsations propres ω_i et des modes propres associés φ_i.

Dans le cas de l’analyse temporelle, on calcule à chaque instant t, par intégration directe des équations d’équilibre, les déplacements aux nœuds q(t) ainsi que les vitesses q'(t) et les accélérations q¨(t).

La seconde approche présente l’avantage de permettre de traiter des sollicitations non stationnaires.

B.1 Analyses reposant sur une recherche modale

B.1 Analyses reposant sur une recherche modale

B.2 Analyses reposant sur une intégration temporelle directe

B.2 Analyses reposant sur une intégration temporelle directe

B.3 Prise en compte de l'amortissement

B.3 Prise en compte de l'amortissement

B.4 Spécificités de l'analyse sismique

B.4 Spécificités de l'analyse sismique

$translationBooks

B.1 Analyses reposant sur une recherche modale

B.1 Analyses reposant sur une recherche modale


Rappel sur la notion d’oscillateur simple – Notion de spectre de réponse élastique

Dans la suite de cette section, beaucoup de méthodes employées renvoient vers la réponse d’oscillateur amorti à un degré de liberté souvent dénommé oscillateur simple.

On rappelle très brièvement les grandes lignes de son comportement ici.

On se propose d’expliciter son comportement en reprenant l’équation de la dynamique appliquée à un oscillateur à un degré de liberté soumis à un chargement temporel harmonique et en bâtissant sa fonction de transfert paramétrée par le ratio entre sa pulsation propre et celle du chargement harmonique ainsi que son taux d’amortissement critique.

On rappelle l’expression de l’équilibre dynamique appliquée au système masse, ressort et amortisseur à un seul degré de liberté.

Qui s’écrit également sous forme canonique en divisant tous les termes par la masse m telle que :

Si on applique un chargement harmonique , on peut trouver la fonction de transfert en réponse absolue en établissant le rapport entre ce chargement d’entrée et la réponse absolue de sortie ci-dessous :

La norme de cette fonction est l’amplitude de la fonction de transfert qui permet de mettre en évidence le phénomène d’amplification dynamique illustrée pour différentes valeurs de taux d’amortissement critique sur la figure 1. L’argument de cette fonction de transfert est la phase.

Fonction de transfert

Un spectre de réponse est différent d’une fonction de transfert. On appelle spectre de réponse élastique la courbe donnant l’accélération (dénommée spectrale) en fonction de la période (ou de la fréquence). Le spectre correspond à l’accélération absolue maximale vue par un oscillateur simple au cours du temps en fonction de sa période propre (ou de sa fréquence propre) et de son taux amortissement critique. Il dimensionne le mouvement sismique. Il est possible de construire une relation approchée entre le spectre de Fourier d’une accélération et le spectre en vitesse spectrale pour un taux d’amortissement nul.

A partir de l’équation de la dynamique de l’oscillateur simple soumis à une accélération quelconque :

On trouve les réponses maximales à partir de sa résolution en temps par une méthode de son choix (intégration de Duhamel, calcul direct …) puis par le calcul des valeurs dites spectrales de déplacement relatif, vitesse relative et accélération absolue :

La notion de pseudo-déplacement, vitesse ou accélération est souvent employée. Il s’agit d’une approximation de ces quantités à partir de l’une d’elle reposant sur le fait que le taux d’amortissement est assez faible (moins de 20%).

Un spectre peut se présenter sous la forme de la figure 2.

Spectre

Notion de base modale

On considère que la réponse de la structure s’appuie sur une combinaison de réponses harmoniques. A ce titre, les réponses harmoniques sont des champs solutions q(t) de l’équation :

Ils peuvent donc s’écrire sous la forme :

Si l’on fait l’hypothèse de Basile (on néglige C), on arrive à rechercher les solutions de l'équation algébrique :

Si l’on ne fait pas l’hypothèse de Basile, cela revient à rechercher les modes complexes ; on postule :

et on résout le problème aux valeurs propres :

Pour y parvenir, sans négliger C, il a fallu poser :

Avec les matrices A et B telles que : 

 et 

On désigne aussi par γ le vecteur tel que :

On recherche les solutions sous la forme de paires d’harmoniques complexes conjuguées telles que :

avec :

  • ψi qui est un vecteur déformée modale complexe,
  • λi qui est une fréquence propre complexe.


Dans tous les cas, réel ou complexe, on a donc à résoudre une équation en λi2 de même degré N (2N pour les cas complexes) que l’ordre de la matrice N (2N pour les cas complexes). L’ordre de cette matrice est égal au nombre de degrés de libertés du système discrétisé (le double pour les cas complexes). L’équation évoquée n’est autre que les cas de nullité du déterminant de cette matrice et donc de son polynôme caractéristique.

A des fins de simplification, on considère par la suite le cas de modes réel avec l’hypothèse de Basile.

A chaque pulsation ωi est associé un vecteur propre. Celui-ci va être recherché en fixant une de ses composantes à 1, puis un système à n-1 paramètres est résolu.

 

Attention : les modes sont donc définis à une constante multiplicative près, d’autant qu’ils sont par la suite normés au moyen de procédés variables. Le plus courant pour les codes EF est de normer les modes au moyen de la matrice masse ; ce point est détaillé par la suite. Ils peuvent être également normés par leur plus grand déplacement modal.


Les modes propres ont la propriété d’être orthogonaux pour la matrice de masse, ce qui se formule comme suit :

En l’absence de chargement, les modes n’ont aucun sens physique. On peut donc choisir de les normer de différentes façons pour rendre leur visualisation compréhensible. A ce titre, il est le plus souvent rencontré dans les codes EF :

  • de normer les modes sur la matrice masse, on a alors :

  • de normer les modes à partir d’un mode en particulier.

Pour les cas où le code ne norme pas cette valeur, on désigne cette grandeur par le terme de masse généralisée et on écrit pour le mode i :

L’objet de ce document étant de traiter des études conduites par les codes EF, on retiendra que les modes sont normés sur la matrice masse et donc la masse généralisée est toujours égale à l’unité.

On construit également la grandeur appelée masse modale qui permet d’identifier la quantité de masse de la structure entraînée par un mode dans une direction donnée Δk:

Attention : la masse modale est une grandeur bien différente de la masse généralisée comme l’équation ci-dessus permet de le mettre en évidence.


Utilisation de la base modale

Réponse temporelle par projection sur base modale

Une solution à l’équation globale :

 

peut être décomposée sur la base de vecteurs propres orthogonaux en N problèmes indépendants qui sont ceux décrivant la réponse d’un oscillateur simple ri(t) :

Pour les cas de séisme, le chargement f(t) est un chargement inertiel appliqué à l’ensemble de la structure : 

.

Analyse harmonique

On s’intéresse dans ce type d’analyse à construire une fonction de réponse en fréquence de la structure aux différents nœuds du modèle. Cette réponse peut être une quantité d’intérêt quelconque (déplacement, vitesse, accélération – absolus ou relatifs –, efforts …) en fonction d’une gamme de chargements harmoniques d’entrée qui peuvent également être décrits de différentes façons semblables à la quantité d’intérêt extraite.

Le chargement imposé sur un ensemble quelconque de nœuds de la structure est dans ce cas précis une harmonique du type : 

.

Le résultat de l’analyse en tout nœud de la structure est une fonction de transfert complexe dont la norme (l’amplitude) et l’argument (la phase) sont le plus généralement utilisés. Ils permettent d’identifier les résonnances de la structure ou de l’équipement et de connaître l’amplitude de sa réponse à différents chargements harmoniques.

Troncature de base modale – Cas général

En règle générale, comme on ne peut pas extraire par calcul tous les modes (trop coûteux et long, voire numériquement inatteignable dans certains cas), il faut se contenter de ceux susceptibles de répondre au chargement de la structure.

Si la fréquence représentative retenue pour le chargement est fc (Hz) on cherchera à extraire les modes jusqu’à 2 fc.

Il est important également de veiller à bien représenter les modes susceptibles de contribuer localement.

Une recherche de mode ne doit pas s’accompagner d’une lecture de valeurs de fréquences mais d’une visualisation de leurs déformées.

$translationBooks

B.2 Analyses reposant sur une intégration temporelle directe

B.2 Analyses reposant sur une intégration temporelle directe

Les schémas d'intégration

Grands principes

Le principe des différentes méthodes d’intégration directe est de découper l’intervalle d’étude en n intervalles de longueur ∆t = T/n et de vérifier l’équilibre aux instants discrets Ti = i ∆t = i T/n. La différence entre les différentes méthodes (différences centrées, Wilson, Newmark…) repose sur l’hypothèse qui est faite sur la variation des grandeurs cinématiques sur l’intervalle ∆t.

Schéma explicite :

Si la valeur de déplacement peut être calculée directement à partir de la ou des valeurs au pas de temps précédents, la résolution est dite explicite. Dans ce cas, l’équilibre est considéré en début d’incrément.

Schéma implicite :

Les méthodes implicites sont celles pour lesquelles l’équilibre doit être considéré à la fin de l’intervalle, ce qui nécessite la résolution d’un système linéaire.

Comme indiqué plus haut, il existe un grand nombre de schémas d’intégration et l’objectif n’est pas ici d’en faire une description exhaustive (voir pour cela des ouvrages spécialisés). Nous nous contentons d’en présenter deux, la méthode des différences centrées, méthode explicite conditionnellement stable et la méthode de Newmark, schéma inconditionnellement stable (pour les problèmes linéaires).

Un schéma est dit inconditionnellement stable si, pour n’importe quelles conditions initiales, la solution reste bornée quel que soit le pas ∆t choisi , en particulier lorsque ∆t/T est grand. Par opposition, une schéma est conditionnellement stable si la solution obtenue reste bornée seulement si ∆t reste inférieur à une valeur limite ∆tcrit.

La précision est un concept différent de la stabilité, qui prend toute son importance pour les schémas inconditionnellement stables. Au-delà des inévitables erreurs d’arrondi, la précision agit sur deux sources d’approximations : un allongement (artificiel) de la période et une diminution de l’amplitude. L’influence de ces deux phénomènes augmente avec l’augmentation de ∆t mais de manière indépendante.

Méthode des différences centrées

Il s’agit d’un schéma de type différences finies reposant sur l’approximation de l’accélération (développement de Taylor au second ordre) :

Pour obtenir le même ordre d’erreur sur la vitesse, on utilise :

Les déplacements à l’instant t+∆t sont obtenus en considérant l’équilibre à l’instant t :

Soit, en introduisant les approximations de l’accélération et de la vitesse :

De la forme :

Cette méthode nécessite une procédure de démarrage afin de calculer q(-∆t) (à partir de l’équilibre à t=0). L’amortissement introduit par ce schéma est nul (pas de décroissance d’amplitude).

Méthode de Newmark

Ce schéma repose sur les approximations suivantes de la vitesse et du déplacement en fin d’intervalle :

De la valeur des 2 paramètres α et δ (compris entre 0 et 1) dépendent la précision et la stabilité de la méthode : le couple (δ = 1/2 et α = 1/4) conduit au schéma inconditionnellement stable appelé méthode de Newmark : il correspond à considérer l’accélération moyenne constante. En considérant l’équilibre à la fin de l’intervalle d’étude (à l’instant t+Δt), on obtient :

avec

et

les coefficients étant définis ci-dessous :

Ce schéma nécessite également une procédure de démarrage : la valeur de q ̈(0) est obtenue en considérant l’équilibre à t=0. Comme la méthode des différences centrées, le schéma de Newmark de base (δ = 1/2 et α = 1/4) n’introduit pas d’amortissement numérique. Le couple (α = 0, δ = ½) permet de retrouver la méthode des différences centrées.

Choix de la discrétisation spatiale et temporelle

Critère sur les tailles d’éléments satisfaisant les longueurs d’ondes

Différents critères peuvent intervenir qui sont directement liés à la précision des résultats attendus et au type de calcul qui est entrepris.

Pour les analyses transitoires, la recommandation généralement retenu est d’avoir de l’ordre de 8 à 10 éléments par longueur d’onde. Les ondes stationnaires étant constitués de la somme d’ondes propagatives, les analyses en dynamique stationnaire (cas des séismes typiquement) verront s’appliquer le même type de critère.

On rappelle l’expression générale reliant longueur d’onde λ, fréquence f et célérité c de l’onde ou encore la pulsation ω est :

En fonction des types d’éléments utilisés et du type d’onde qui retient notre intérêt, on retiendra les différentes formulations suivantes :

Éléments volumiques dans un milieu isotrope (cas type d’une modélisation d’un sol élastique de module E, de coefficient de Poisson ν et de masse volumique ρ) :

  • Onde de cisaillement : 
  • Onde de compression : 

Éléments coques isotropes :

  • Onde de flexion :  où h est l’épaisseur de la coque, D le coefficient de raideur flexionnel.
  • Onde de cisaillement : 
  • Onde de compression : 

Ces formulations s’écrivent plus généralement de la façon suivante pour un milieu anisotrope avec 2 directions principales 1 et 2 :

  • Onde de flexion : 
  • Onde de cisaillement :  (expression valable en coque et plaque et en volume),
  • Onde de compression :  (expression valable en coque et plaque et en volume).


On peut trouver la taille requise pour cibler une fréquence de coupure de 40 Hz dans des éléments coques et plaques ou dans des éléments volumiques. La taille requise est donnée par le nombre ci-dessous / λ pour 10 éléments.

Types de modélisations Fréquence cible [Hz] Epaisseur [m] Type d'onde Coeff de taille d'élément
Coques et briques - GC Béton 40 0,2 Flexion 0,59
40 0,25 Flexion 0,66
40 0,5 Flexion 0,93
40 1 Flexion 1,32
40 1,5 Flexion 1,61
40 - Cisaillement 6,04
40 - Compression 6,75
Coques et briques - Méca Acier 40 0,01 Flexion 0,16
40 0,02 Flexion 0,22
40 0,1 Flexion 0,50
40 - Cisaillement 8,02
40 - Compression 9,58


Ce critère n’est pas nécessairement suffisant et, pour les cas où une analyse temporelle par projection sur base modale ou une réponse spectrale est recherchée, une analyse de sensibilité doit être mené afin d’observer si un raffinement de maille conduit à faire varier significativement les facteurs de participation modaux (ou alors la masse modale si ce critère est scruté). Ces éléments sont décrits plus loin dans ce chapitre et leur lien sont exposés.

Critère sur les pas de temps

Les schémas temporels conditionnellement stables comme le schéma des différences centrées doivent satisfaire une condition sur le pas de temps choisi. On appelle souvent cette condition la condition CFL (Courant-Friedrichs-Levy).

Cas non amorti Cas amorti
Satisfait pour  Satisfait pour 

Exemple d’application pour un modèle de treillis métallique dont le plus petit élément est de longueur L :

On considère un élément barre à 2 nœuds notés 1 et 2.

Les barres sont de longueur L, de section A et de masse volumique ρ et de module d’Young E.

On considère la matrice masse telle que la masse totale est répartie également aux deux extrémités.

Les matrices K et M s’expriment respectivement sous la forme :

 et 

Le polynôme caractéristique s’exprime tel que :

Ce qui conduit à :

Comme la vitesse de propagation des ondes de compressions dans un milieu continu s’exprime par :

On trouve que :

 et  soit 

Dans un modèle dont le maillage n’est pas régulier, c’est le plus petit pas de temps qui pilote le pas temps global.


$translationBooks

B.3 Prise en compte de l'amortissement

B.3 Prise en compte de l'amortissement

L’amortissement des amplitudes des phénomènes vibratoires d’une structure est lié à des phénomènes dissipatifs d’origines très variées qui peuvent se produire au sein de la structure. Il peut s’agir de phénomènes visqueux propres à des amortisseurs fluides que l’on peut trouver dans des boîtes à ressort ou des amortisseurs (isolateurs à câble, amortisseur de type JarretTM, amortisseur à fluide visqueux). Il peut s‘agir également de phénomènes de dissipation structurelle liés à des états de fissurations et aux frottements aux interfaces des matériaux. L’énergie est consommée dans les hystérésis de comportement des matériaux.

Dans l’essentiel des cas, la prise en compte de l’amortissement s’opère dans les calculs dynamiques au travers d’une matrice d’amortissement visqueux C. Les termes de cette matrice peuvent être localement définis exactement pour les cas où l’amortissement est totalement maîtrisé parce qu’il s’agit d’équipements amortissants dont on connaît les caractéristiques. Mais la plupart du temps, pour les structures de GC, il est rare de connaître le terme d’amortissement compte tenu de l’hétérogénéité des structures et du caractère non visqueux des matériaux eux-mêmes.

On préfère construire cette matrice à partir des taux d’amortissement critique, définis réglementairement pour les différents matériaux. Ces taux sont définis avec une hypothèse de comportement visqueux proportionnel à la vitesse de sollicitations ; il s’agit d’une fraction de la valeur d’amortissement critique qui conduit pour un oscillateur simple à un retour à l’équilibre sans passer par un régime d’oscillations.

Pour les cas où les termes d’amortissement sont connus, la matrice C est construite explicitement et il est nécessaire d’avoir recours à une méthode de calcul temporel dédiée, afin de bien respecter les caractéristiques de celle-ci :

  • soit avec un calcul par intégration directe : il s’agit de résoudre l’équilibre dynamique à chaque instant en ayant discrétisé le temps,
  • soit avec un calcul par projection sur base modale, faisant appel à des modes complexes.

Les analyses de phénomènes de dynamique rapide (explosion, impact) doivent faire l’objet de considérations spécifiques.

En effet, il n’est pas nécessaire de considérer d’amortissement pour les études de vérifications locales de résistance ou de vibrations sur des temps trop courts pour « activer » des réponses harmoniques de la structure. On peut donc se passer pour ces cas de prendre en compte l’amortissement.

Il convient par ailleurs de noter que ce genre de calculs, effectués par des codes de dynamique faisant appel à des schémas d’intégration temporelle dit explicite (LS-Dyna, Radioss, PAM-Crash, Europlexus, …) n’ont pas recours à une résolution classique faisant intervenir une « inversion » de la matrice de rigidité. Pour les analyses de vibrations induites sur une fenêtre de temps succédant la durée du choc, il est par contre nécessaire de caractériser l’amortissement.

Nota : dans des calculs de dynamiques rapides (quelques millisecondes), le terme d’amortissement est généralement négligé car son ordre de grandeur est très inférieur à celui des termes inertiels et de rigidité locale sur le temps très court où le phénomène est analysé.


Pour les cas plus « classiques », on a généralement recours à une construction de la matrice d’amortissement à partir des taux d’amortissement critique définis pour les matériaux et par combinaison linéaire des matrices M et K qui ont déjà été élaborées par le code de calcul aux éléments finis. Cette démarche permet de simplifier considérablement les étapes de calcul. Cette approche est liée à l’hypothèse dite « de Basile » pour les analyses modales. Cette hypothèse a été formulée dans le cadre de l’identification de l’amortissement structurel (principalement en aéronautique), et est restitué comme suit :


« Même en présence de couplage modal d’amortissement, les équations modales du mouvement sont dynamiquement découplées, pour les structures faiblement amorties, si la séparation des modes en fréquences est satisfaisante ».

Attention : prendre en compte l’hypothèse de Basile conduit à sous-évaluer la réponse dynamique ailleurs qu’à proximité des matériaux ou points d’appuis qui présentent un amortissement élevé par rapport à un calcul temporel par intégration directe ou à l’utilisation de méthodes plus sophistiquées comme celles des modes complexes. Ceci justifie le fait d’avoir recours à un calcul par intégration directe des équations de mouvement ou en faisant appel à des modes complexes.


Le plus souvent, l’hypothèse de Basile est considérée et la matrice est construite à partir de la méthode de Rayleigh.

Les termes α et β sont à caler. De manière pragmatique, la démarche qui suit doit être prise en compte. On considère :

, soit exprimé en fréquence : ,

, soit exprimé en fréquence : ,

f1 étant une valeur inférieure au premier mode significatif de la structure, et fn le premier mode rencontré après la coupure du spectre de réponse élastique caractérisant le signal temporel appliqué. On reconstitue pour toute valeur de ω la valeur de l’amortissement réduit, en faisant appel à la formule ci-dessous :

, soit exprimé en fréquence : 

Attention : l’amortissement étant très largement surévalué au-delà des 2 fréquences d’intérêt, il est nécessaire de maîtriser ces deux valeurs. Si fn a bien été choisie, la surévaluation de l’amortissement est sans conséquence car on trouve sur ces plages de fréquences des modes de corps rigides qui seront sollicités dans la réponse et les modes au-delà de la fréquence de coupure sont insensibles à l’amortissement. Mais il faut être sûr pour cela que l’on ne prête pas attention à d’éventuels modes locaux à haute fréquence.

L’impact d’une erreur dans le choix de la première fréquence peut être plus lourd de conséquences si un mode significatif a été oublié car sa réponse sera dès lors négligée.

Il importe de noter qu’entre les deux fréquences, l’amortissement est sous-évalué. On a en ce cas une évaluation conservative de la réponse de la structure.



$translationBooks

B.4 Spécificités de l'analyse sismique

B.4 Spécificités de l'analyse sismique


Réponse spectrale - Cas spécifique du séisme

Le principe de la méthode est, pour une direction de séisme de donnée, de construire les réponses maximales à partir du spectre de chargement en tout point considéré mode par mode, puis de les cumuler entre elles par différentes méthodes.

Une fois obtenue la réponse sismique pour une direction donnée, les directions sismiques sont combinées entre elles pour avoir la réponse globale.

On se place dans un cas où les vecteurs modaux sont normés au moyen de la matrice de masse, i.e. dans le cas de la plupart des codes de calcul aux éléments finis.

On peut écrire l’expression pour un degré de liberté sous la forme canonique :

Le terme pik est appelé coefficient participatif. Il s’agit d’une notion importante dans le cas de l’analyse sismique car elle restitue la contribution physique d’une excitation dans une direction k donnée.

Elle est déterminée pour les cas où l’on norme les vecteurs propres sur la matrice masse par :

D’autre part, on connaît pour un mouvement sismique  donné, dans une direction sismique , les solutions majorant en moyenne les valeurs maximales de  et  en termes de pseudo-accélérations, pseudo-vitesses et déplacements via les spectres respectifs de ce mouvement :

Cette réponse est construite pour un amortissement modal ξi qui peut être évalué au prorata de l’énergie de déformation modale (voir §1.2.1 de Génie parasismique Tome 3, J. Betbeder-Matibet).

où :

  • EiT est l’énergie de déformation totale de la structure pour le mode i telle que : 
  • Ei matériau j est l’énergie de déformation dans la structure dont l’amortissement matériau est ξj.

Ensuite à partir de la réponse de l’oscillateur simple, on reconstruit la réponse pour chaque mode, respectivement en accélération, vitesse et déplacement, dans la structure au moyen de :

Ces grandeurs modales doivent ensuite être combinées entre elles pour une direction de chargement donnée afin d’obtenir la réponse complète de la structure. Il existe diverses méthodes de recombinaison (par la suite la notation ak,p désigne à titre d’exemple la composante p du vecteur ai,k, il en va de même pour φi,p à partir de φi) :

  • Racine carrée de la somme des carrées, dit « cumul quadratique simple » ou également « SRSS : Square Root of the Sum of the Square » – En chaque nœud p :

  • Cumul quadratique complet (dit CQC) – En chaque nœud p :

où le terme Pij désigne le coefficient de couplage quadratique des modes i et j. Il s’écrit :


Il importe de noter que, dans cette combinaison, il s’agit d’un cumul de termes algébriques mais dont le signe est toujours positif et ne pose donc pas de difficulté à être placé sous la racine.

En considérant des amortissements modaux ξi et ξj égaux à une valeur variable telle que x = {0.2, 4, 7, 20 et 30%}, on peut tracer pour des ratios de pulsations ωij le graphique de la figure 3 qui permet de faire apparaître que les modes voisins se cumulent très nettement tandis que le coefficient de cumul diminue assez rapidement lorsque ce ratio croît et ceci d’autant plus lorsque l’amortissement est faible.

Illustration du coefficient de couplage pour des taux d’amortissements variables en fonction des ratios de pulsation

Ce dernier constat milite pour l’utilisation d’un cumul CQC lorsque des modes sont assez proches car le cumul SRSS ne prendrait pas en compte le cumul des modes proches.

En revanche, il n’y a aucune plus-value à utiliser un cumul CQC par rapport à un SRSS lorsque les modes sont éloignés les uns des autres.

 

Attention : Les réponses modales cumulées ne doivent pas être utilisées pour calculer d’autres quantités d’intérêt. A titre d’exemple il ne faut pas surtout pas évaluer les efforts issus des réponses CQC en calculant  : cette évaluation est erronée par rapport à l’évaluation par cumul des efforts modaux.

Il faut donc pour cet exemple précis dérouler le calcul comme suit :

  1. ,

  2. Pour chaque terme p du vecteur fk, on écrit : ,

Il est aussi important de garder à l’esprit que de la même manière des critères permettant d’évaluer des états de fissurations sur la base des invariants du tenseur des contraintes, et donc des contraintes principales, (cf. par exemple le critère proposé dans l’annexe LL de l’EN1992-2) ne peuvent être considérés au terme d’un calcul avec un cumul quadratique de quelque nature que ce soit.

Les réponses spectrales de chaque mode une fois cumulées pour une direction de chargement sismique doivent être combinées afin d’obtenir la réponse totale. On parle dans ce cas de cumul spatial.

Les cumuls spatiaux peuvent s’opérer selon diverses méthodes :

  • Par racine carré de la somme des carrés des réponses obtenue dans chaque direction. Cette méthode fait perdre le signe et toute concomittance logique des sollicitations. Elle fournit une seule quantité scalaire pour chaque quantité d’intérêt :

  • Par cumul algébrique de type « Newmark ». Cette approche repose sur une hypothèse d’indépendance les unes aux autres de chaque réponse spatiale et pondère d’un coefficient μ, dont la valeur varie suivant les standards, les 2 autres réponses défavorablement par rapport à celle d’une direction privilégiée. Elle prend en compte la variabilité de signe de chaque quantité d’intérêt. Elle conduit donc non pas à une seule réponse cumulée mais à 24 au total comme cela est compréhensible dans l’équation suivante :

Troncature de base modale – Cas du séisme

Pour le cas spécifique du séisme et en retenant une norme des modes par la matrice masse, il vient :

Il s’agit de la masse entraînée pour un mode i dans une direction k. Cette masse est donc liée à une direction spécifique, pour un mode i il y aura donc 3 masses modales en fonction des différentes directions de sollicitation (si on est dans l’espace, 2 dans le plan et une seule pour un problème unidimensionnel), la figure 4 apporte une illustration en 2D pour un portique.

Les déformées modales multipliées par les coefficients participatifs représentent la déformée de la structure dont le produit en tout point avec la réponse d’un oscillateur simple fournit la réponse au cours du temps pour un mode déterminé. Ceci est détaillé dans le chapitre suivant.


Illustration schématique des entraînements de masse pour un même mode suivant 2 directions différentes

Comme on l’a évoqué plus haut si le modèle possède N degrés de libertés, il y aura au plus N modes propres.


Toutefois, dans la mesure où la recherche des fréquences et modes propres fait l’objet d’une démarche numérique, tous les modes ne sont pas extraits. En théorie, on doit avoir :

, ou  ou 


En pratique, le code va s’arrêter à une limite fixée par l’utilisateur au travers d’une fréquence maximale d’extraction. Puis l’utilisateur doit vérifier qu’il en dispose d’un nombre suffisant pour restituer un pourcentage de masse fixé par la règle d’ingénierie qu’il applique. Le plus généralement ce critère est de 90 % de la masse.

Attention : Si, à la fréquence de coupure, le pourcentage de masse ciblée n’est pas atteint, il peut être nécessaire de faire appel à un pseudo-mode.

Inversement si un pourcentage de masse important est atteint à une fréquence très basse, très en-dessous de la fréquence de coupure, il est nécessaire de faire appel à un pseudo-mode ou de compléter la base modale. En effet, la masse participante d’un plancher pour un mode de flexion local sur un très grand ouvrage représente une très faible fraction de la masse totale.

Attention également aux modes antisymétriques dont la masse modale peut être nulle car les masses entraînées autour d’un axe de la structure se compensent.

Pseudo-mode ou correction statique

Comme on l’a formulé précédemment, le critère de sélection des modes propres porte sur le cumul de la masse modale fixé à 90 % pour une fréquence de l’ordre de la fréquence de coupure du spectre sismique i.e. 40 Hz max.

Lorsque cette valeur n’est pas atteignable, le pourcentage plus faible retenu va être complété d’une masse additionnelle associé à un « pseudo-mode » de vibration.

On rappelle que l’on a pour une direction de séisme k une réponse cumulée sur n modes :

En toute rigueur, s’il y a N degrés de libertés, cette formulation serait :

Comme l’on a vu précédemment, les modes au-delà de la fréquence de coupure sont des modes rigides. La structure réagit en phase avec le chargement sismique qui lui est appliqué avec un déplacement relatif nul.

Pour compléter la base modale, on construit les pseudo-modes en considérant que la réponse totale est la somme de la réponse « dynamique » prenant en compte la base modale sur les n premiers modes retenus et d’un terme proportionnel à l’accélération sismique du support q̈s(t). On retient les expressions suivantes :

où Pk est la déformée dans la direction k de la structure soumise au chargement statique équivalent à la masse de la structure accélérée à l’accélération correspondante au dernier mode extrait de pulsation n :

Comme la réponse en accélération relative pour la structure est nulle au-delà de la fréquence de coupure, on peut simplement écrire :

En reprenant la formule précédente en fonction de , on pourra plus facilement faire le lien dans le cadre des analyses spectrales avec le spectre de pseudo-accélérations :

On considère donc également la formule :

Une fois ces grandeurs vectorielles évaluées, se pose la question de leur cumul.

En effet, la solution temporelle pour une direction de séisme est donnée via la recombinaison des composantes sur la base modale :

En revanche, cette combinaison n’est pas applicable sur les grandeurs maximales que nous venons d’évaluer.

$translationBooks

Chapitre C - Calculs statiques non-linéaires

Chapitre C - Calculs statiques non-linéaires

Les calculs non-linéaires sont généralement longs : un programme d’éléments finis les résout comme une succession de problèmes linéaires enchaînés (processus incrémental). La préparation de tels calculs demande à l’utilisateur de faire des choix, comme la définition des incréments, le choix d’un algorithme, etc. Le traitement exige donc une certaine expérience d’utilisation du logiciel.

Avant de se lancer dans un calcul non-linéaire, il importe au préalable de préciser ce que l’on en attend : cela permet en effet de savoir quand arrêter le calcul ; car la plupart des calculs non-linéaires peuvent être a priori poursuivis jusqu’à l’apparition d’un ou plusieurs mécanismes ou, dans le cas de problèmes de contact, jusqu’à la disparition du nombre minimum de liaisons (ou d’appuis).

C.1 Les problèmes de mécanique non-linéaire

C.1 Les problèmes de mécanique non-linéaire

C.2 Pourquoi exécuter des calculs non-linéaires ?

C.2 Pourquoi exécuter des calculs non-linéaires ?

C.3 Mise en œuvre

C.3 Mise en oeuvre

C.4 Problèmes de convergence ? Symptômes et solutions

C.4 Problèmes de convergence ? Symptômes et solutions

$translationBooks

C.1 Les problèmes de mécanique non-linéaires

C.1 Les problèmes de mécanique non-linéaires


1) Description des non-linéarités potentielles

La non-linéarité d’un problème de mécanique vient de ce que les coefficients apparaissant dans l’équation d’équilibre dépendent de la déformation du solide à l’équilibre. Autrement dit, l’équation d’équilibre est généralement une équation de type implicite.

Il existe plusieurs catégories de non-linéarités pour les problèmes de mécanique statique :

  • Les non-linéarités matérielles : c’est le cas dès lors que la loi de comportement n’est pas linéaire ou que la réponse dépend de l’histoire du chargement. Autrement dit, la contrainte n’est pas une fonction linéaire de la déformation. L’exemple le plus fréquent en génie civil est celui des matériaux sollicités au-delà de leur limite élastique, et qui développent un comportement élasto-plastique. Ce comportement se caractérise par une dépendance de la raideur du matériau vis-à-vis de son état de contrainte.
  • Les non-linéarités géométriques : c’est le cas dès lors qu’on travaille en grands déplacements ou en grandes déformations. Dans le premier cas, on ne peut plus écrire le problème en négligeant les changements de configurations. Dans le second cas, on ne peut plus approcher la déformation par un simple gradient du déplacement.
  • Les non-linéarités des conditions aux limites : cela peut être le cas, lorsque chargement est progressif, lorsqu’il y a contact potentiel entre 2 corps, dans le cas de force suiveuse, mais aussi lorsqu’on simule le phasage de construction, ou le lançage d’un tablier, le creusement d’une galerie, la mise en remblai, etc.

Ces types de non-linéarités peuvent être couplés, si le code de calcul le permet, mais cela complexifie d’autant la résolution du problème.

2) Principe de résolution d'un problème non-linéaire : Méthode de Newton

Lorsqu’on résout le problème élément fini, on cherche le champ de déplacement u, tel que les efforts intérieurs Lint soient égaux aux efforts extérieurs Lext :

, qui est un problème non-linéaire en fonction de u.

Pour résoudre un problème de statique non-linéaire, on utilise généralement un algorithme incrémental. Pour cela le problème est paramétré en fonction d’un paramètre t (qui est un pseudo-temps, contrairement aux calculs dynamiques). Ce paramètre sert à indexer les chargements successifs imposés à la structure. Plus précisément, on va chercher l’état d’équilibre correspondant aux chargements successifs F1, F2, …

Ce découpage conduit à résoudre une succession de problèmes quasi-linéaires comme illustrée sur la figure ci-dessous, et à déterminer l’état de la structure à l’instant t (déplacements, déformations, contraintes) en connaissant la solution à l’état t-1. Plus le nombre de découpage est important, meilleure est la précision.

Principe du paramétrage de problème en fonction de t

A chaque incrément ti, le problème discrétisé à résoudre est : Ki x qi =Fi où qi est le vecteur déplacement inconnu sous le chargement imposé Fi. Alors que dans le cas linéaire vu au chapitre 1, la matrice K était explicite, lorsque le problème est non-linéaire, Ki est une matrice dont les termes dépendent implicitement de la valeur de qi. Déterminer qi ne peut donc pas se faire directement en inversant la matrice K.


La méthode la plus courante pour résoudre cette équation non-linéaire est d’utiliser un algorithme de type Newton. Le principe est de construire une bonne approximation de la solution de l’équation  en considérant son développement de Taylor au premier ordre . On part donc d’un point initial (suffisamment proche de la solution) et on construit par récurrence la suite .


A chaque itération, on évalue le vecteur résidu F(qk) et on poursuit tant que celui-ci n’est pas inférieur (en norme) à une valeur arbitrairement proche de 0. Ce critère de convergence doit être choisi avec soin en fonction de la norme utilisée par le code de calcul (voir §3.3 pour plus de détails).


Remarque :
 Avec la méthode de Newton, à chaque itération, on doit calculer la matrice tangente au point considéré .

Le coût de calcul de cette matrice peut être coûteux en temps. Si l’utilisation de cette matrice permet d’avoir une convergence quadratique (donc en peu d’itérations), il n’est pas indispensable d’utiliser cette matrice. D’autres stratégies peuvent être appliquées pour estimer cette matrice. On parle alors de méthode de quasi-Newton. Ainsi il est ainsi envisageable d’utiliser la matrice tangente sans la réactualiser à chaque itération mais aussi d’utiliser la matrice élastique (figure b) ou la matrice sécante dans les cas de modèle endommagement. Une illustration des itérations successives selon la matrice utilisée est proposée sur la figure ci-dessous.

Illustration de la méthode de Newton ou quasi-Newton (matrice élastique)

En règle générale, l’utilisation de la matrice tangente permet une convergence plus rapide (en moins d’itérations), mais les variantes peuvent être plus efficaces ou plus robustes selon les situations.

La méthode de résolution étant itérative, le processus doit être interrompu lorsque le critère d’arrêt est atteint, c’est-à-dire que l’on a vérifié qu’une grandeur (plusieurs définitions sont possibles) devient négligeable. L’algorithme global est le suivant :

en désignant les incréments, i indexant les itérations de Newton et ε étant une valeur positive, proche de 0.

Remarque : l’algorithme de Newton est utilisé pour résoudre l’équilibre à chaque pas de temps. Il peut être également utilisé pour trouver la contrainte en chaque point de Gauss (à toutes les itérations du problème de Newton à l’échelle globale) lorsque la loi de comportement le nécessite.

$translationBooks

C.2 Pourquoi exécuter des calculs non-linéaires ?

C.2 Pourquoi exécuter des calculs non-linéaires ?

Comme dans bien des domaines de la physique, il n’est raisonnable d’entreprendre un calcul non-linéaire qu’une fois que l’on a au moins une petite idée de la « fin de l’histoire », c’est-à-dire de la façon dont la structure va évoluer au point de devenir instable. Nous présentons ci-dessous les raisons pouvant justifier un calcul non-linéaire, en distinguant les bonnes de certaines plus discutables…

1) Bonnes raisons
  • Étudier la redistribution des efforts. Une fois que certaines parties d’une structure sont entrées en plasticité, le niveau de contrainte y est, en un certain sens, à peu près « gelé » à une valeur déterminée, et le travail des actions extérieures ne peut plus être dissipé qu’en augmentant le niveau de contrainte ailleurs, ou en faisant s’écouler la matière des zones déjà plastiques. Le but est ici de vérifier dans quel ordre les éléments d’une structure « lâchent », et d’appréhender le mécanisme de ruine final.
  • Déterminer l’évolution des conditions d’appui d’une structure, soit par un calcul de contact, soit par formation de rotules plastiques internes.
  • Sans aller jusqu’à s’intéresser à la redistribution des efforts, obtenir un état d’équilibre où, dans les zones plastiques, les contraintes sont, en moyenne, proches de la contrainte d’écoulement du matériau, et où les contraintes sont bien dans le domaine élastique ailleurs. Pour les géomatériaux, obtenir un état d’équilibre où les contraintes sont, sauf en quelques zones très ponctuelles, des contraintes de compression (éliminer les tractions).

Ces démarches sont raisonnables, car :

2) Raisons discutables
  • connaître l’aspect de la structure très déformée (et même en train « d’exploser »), éventuellement en utilisant l’algorithme de remaillage automatique fourni avec le logiciel,
  • faire un calcul avec un comportement a priori très riche du point de vue descriptif, pour voir si ça sert à quelque chose, ou si on comprend mieux le résultat qu’avec le calcul élastique fait précédemment.

Ces motivations sont de celles qui aboutissent le plus souvent à des problèmes de non-convergence ou de « pivot nul » (cf. ci-après). Lorsque l’on interrompt le calcul, faute de mieux, à une étape du processus itératif, la déformée de la structure (donnée par la déformée du maillage, ou le champ des déformations) est illusoire :

$translationBooks

C.3 Mise en œuvre

C.3 Mise en œuvre

La bonne méthode consiste toujours à procéder par étape et à ne pas introduire toutes les non-linéarités simultanément :

  • faire un premier calcul avec une loi élastique avant d’utiliser une loi non-linéaire ;
  • faire un calcul sans contact avant de le mettre en œuvre, faire un calcul avec frottement nul avant d’inclure le frottement, …

Toutes les étapes doivent être vérifiées proprement et il est nécessaire de vérifier que la solution converge lorsqu’on raffine le maillage et le pas de temps.

Enfin, il faut avoir à l’esprit que, dans les cas non-linéaires, il est possible que le problème à résoudre ne possède pas une solution unique (instabilité, bifurcation, etc.).

1) Choix de la discrétisation

La discrétisation doit être choisie avec soin de manière à être cohérente avec les chargements que l’on souhaite appliquer, notamment si ceux-ci ne sont pas monotones.

Plus les non-linéarités sont importantes, plus il faut utiliser des incréments de chargement petits. De même, si le modèle de comportement est complexe, il arrive que l’intégration de la loi de comportement ne soit pas totalement implicite : dans ce cas, il faut vraiment vérifier que les incréments sont suffisamment petits pour que le modèle ait convergé.

Selon les codes de calcul, une gestion plus ou moins automatique de ces pas de temps peut être proposée qui peut permettre de diminuer le pas de temps en cas de problème de convergence voire de l’augmenter si les non-linéarités sont faibles et que la convergence est très rapide. En toute rigueur, il convient de vérifier que la convergence en temps est bien atteinte en effectuant un deuxième calcul en diminuant les incréments de chargement.

2) Choix des conditions aux limites

Il faut avoir conscience que, dans le cas d’une loi de comportement de type adoucissante ou s’il existe une charge limite, le chargement à force imposée peut devenir illicite, comme illustré sur la figure ci-dessous.

Dans le premier cas, plus la force imposée s’approche de la charge limite, plus la convergence sera difficile, jusqu’à être impossible, si l’on cherche à la dépasser. Et dans le deuxième cas, il sera impossible de dépasser l’effort au pic et d’obtenir la solution post-pic.

 

Exemples de réponse avec charge limite ou perte de rigidité

Il existe également d’autres cas, où il n’existe pas une solution unique à effort donné ou déplacement donné. C’est le cas, par exemple dans le cas de flambement de coque mince ou lorsque le problème est tel qu’il existe des bifurcations, par exemple une branche dissipative et une branche élastique (cf. figure ci-dessous).

 

Réponse (a) avec branches multiples (b) non-monotone en déplacement et chargement

Si la modélisation correcte nécessite d’appliquer une force imposée, le problème peut être résolu via des méthodes de continuation ou de longueur d’arc (cf. documentation Code_Aster [R5.03.80] par exemple).

Pour qu’un code de calcul puisse travailler efficacement, il est donc nécessaire de concevoir son modèle, et de préparer les données, de telle façon qu’on travaille à la recherche d’une solution unique, en introduisant diverses restrictions, comme :

3) Évaluation de la convergence de l'algorithme global non-linéaire

Les critères d’arrêt pour l’algorithme de Newton sont à considérer avec attention. Il convient de vérifier avec soin quel critère d’arrêt est programmé par défaut ou proposé par le logiciel utilisé. Par exemple, certains algorithmes vont tester la convergence sur la norme du résidu (dans ce cas, le critère dépend de l’unité de force utilisée dans le calcul, la précision n’est pas la même si le résidu requis est 1N ou 1MN), d’autres proposeront un critère d’arrêt où la norme du résidu est rapporté à la norme du chargement (dans ce cas, on s’affranchit de l’unité physique choisie). Des précautions toutes particulières doivent être prises dans les cas où on mélange différentes modélisations (massifs, éléments de structures) ou lorsqu’on utilise des éléments mixtes.

De manière générale, il faut faire très attention lorsqu’on dégrade le critère de convergence par rapport à la valeur recommandée pour aboutir à un résultat.

$translationBooks

C.4 Problèmes de convergence ? Symptômes et solutions

C.4 Problèmes de convergence ? Symptômes et solutions

À la différence des calculs d’élasticité linéaire, les calculs non-linéaires peuvent effectuer des itérations apparemment indéfiniment, sans jamais vérifier le critère d’arrêt ; mais ils peuvent aussi s’interrompre brutalement à un incrément de chargement donné (au-delà du premier), avec signalement d’un « pivot nul ».

Examinons les causes de ces types d’échec.

1) Non-convergence des itérations

Ce phénomène est presque toujours une conséquence directe de l’inadéquation de la méthode itérative utilisée : le problème à résoudre ne vérifie certainement pas les conditions mathématiques qui assurent sa convergence (opérateur non-contractant, opérateur non-positif, opérateur non convexe, restrictions insuffisantes de l’intervalle dans lequel on recherche la solution, etc.).

Il peut arriver, il est vrai, que ce phénomène résulte d’une accumulation d’erreurs d’arrondis numériques, imputables à la programmation elle-même du code de calcul ; ce n’est toutefois pas souvent le cas.

L’examen des résidus au fil des itérations peut renseigner sur la nature du problème :

  • les résidus tendent bien vers 0, mais de plus en plus lentement. Ce comportement traduit un maillage inadapté au calcul de l’équilibre : c’est le cas dans les phénomènes de localisation, où les non-linéarités se concentrent dans certaines parties de la structure. Dans ce cas, il vaut mieux reprendre le calcul en réfléchissant à un maillage mieux adapté au mécanisme de la structure.
  • les résidus deviennent de plus en plus grands. C’est l’instabilité caractéristique : il n’existe plus, pour cette structure, d’état d’équilibre à ce niveau de chargement.

Les résidus, à partir d’un certain stade, remontent brutalement à une certaine valeur, puis reprennent leur décroissance, etc. Ce comportement survient lorsque l’algorithme tombe sur plusieurs solutions également viables sans autre précision : il converge vers l’une de ces solutions, puis « rebondit » vers une autre. Autrement dit, sans que l’on en ait conscience, il existe plusieurs équilibres possibles !

2) Arrêt brutal

L’arrêt brutal du programme dans un calcul non-linéaire se signale souvent par un « pivot nul » lors de la résolution d’un système linéaire. Cette expression renvoie au « pivot de Gauss », qui est l’algorithme utilisé (aux variantes près) pour résoudre les systèmes linéaires. On a vu plus haut, en effet, que la plupart des algorithmes itératifs remplacent le problème à résoudre par un système linéaire.

Cette situation correspond à une structure devenue instable : en élastoplasticité ou en endommagement, cela signifie que la structure a perdu toute cohésion, aux arrondis de calcul près ; dans les problèmes de grands déplacements ou de contact, cela traduit l’apparition d’un mécanisme, c’est-à-dire la possibilité d’un mouvement de corps rigide : tous les nœuds d’interface sont glissants, et les appuis fixes sont en nombre insuffisant pour obtenir un équilibre.

Il arrive aussi que le programme s’interrompe pour des raisons n’ayant aucun rapport avec l’état de la structure : débordement d’espace-disque, voire bug de programmation. Le premier cas est simple à écarter.

3) Quand rien ne marche

Si un calcul non-linéaire ne converge pas, et que l’on n’a pas trouvé dans les explications ci-dessus d’indications utiles, il y a lieu de passer attentivement en revue les questions suivantes :

$translationBooks

Chapitre D - Génie civil

Chapitre D - Génie civil


D.1 Les matériaux du génie civil

D.1 Les matériaux du génie civil

D.2 Les catégories d'éléments d'ouvrage

D.2 Les catégories d'éléments d'ouvrage

D.3 Les phases de construction

D.3 Les phases de construction

$translationBooks

D.1 Les matériaux du génie civil

D.1 Les matériaux du génie civil

Ce chapitre a pour objectif de présenter les particularités de ces matériaux dans la perspective d'une modélisation aux éléments finis. Parmi celles-ci, nous trouvons les principaux matériaux utilisés en génie civil (béton, bois, métal) en eux-mêmes mais aussi leur utilisation dans des structures.

1) Béton

Le béton est un matériau composite composé de granulats de diverses tailles et d’une pâte. Cette pâte constituée de ciment, d’eau et le plus souvent d'adjuvants, joue le rôle d’une colle (liant) entre les granulats. Le béton bien qu’hétérogène a su révolutionner le monde de la construction grâce à ses propriétés telles que sa souplesse d’emploi, sa moulabilité offrant une variété d’aspects (forme, teinte, texture), sa résistance mécanique, sa résistance au feu, son isolation acoustique et sa diversité (béton armé, béton précontraint, etc.).

Cependant, le béton est soumis à des processus de dégradation. Ces derniers peuvent être issus de sollicitations environnementales telles que l’humidité, la pluie, le froid, le gel-dégel, la chaleur, le vent, la sécheresse, d’agressions provenant de sels de déverglaçage, d’alcali-réaction, de réactions sulfatiques interne et externe ou/et de corrosion et de sollicitations mécaniques comme des chocs, l’augmentation de la charge. Selon le phénomène de dégradation, la détérioration de l’élément structurel est progressive, une fois la contrainte de traction du béton atteinte. Cette dégradation conduit à la coalescence de microfissures vers la formation de macrofissures et la localisation des déformations. L’augmentation de la déformation entraîne la chute de la capacité résistante du matériau ce qui montre un comportement adoucissant ou de softening qui se termine par la ruine de la structure.

Loi de comportement - Dans les calculs par éléments finis, on ne considère généralement que certains aspects du comportement du béton car, c’est un matériau complexe. Pour que la loi de comportement du béton soit la plus pertinente, il est nécessaire de bien comprendre le comportement et les propriétés du béton. De ce fait, non seulement il est essentiel de connaître non seulement le comportement de la pâte de ciment, mais aussi les granulats, qui donnent au béton ses propriétés de rigidité et de raideur. Elle est responsable des propriétés de résistance mécanique du béton. Grâce à l’eau contenue dans la pâte, elle est également la cause de ses principales défaillances telles que l’augmentation de la porosité, la diminution de la résistance mécanique, les effets différés (retrait, fluage) et le transfert d’agents agressifs.

Le béton est considéré comme un matériau quasi-fragile de par son comportement complexe. Les principales caractéristiques de comportement sont mises en évidence dans la littérature. La caractéristique la plus significative propre au béton est son caractère fragile en traction et plus ductile en compression et sujet aux phénomènes unilatéraux lors de cycles traction/compression. D’autres aspects sont à prendre en compte tels que les déformations permanentes et les effets hystérétiques liés aux frottements entre les lèvres des fissures lors de chargements cycliques. Certaines propriétés sont à considérer également comme le retrait et le fluage selon l’environnement et le chargement appliqué au matériau béton.

Dans la littérature, les deux principales approches pour modéliser le comportement du béton sont :

  • Les modèles basés sur une approche continue. Cela veut dire que le béton est considéré comme un milieu continu (Bazant, 1979). Et pour tenir compte de la fissuration du béton, cette dernière est déduite de relations entre les déformations et les contraintes. Il existe les modèles élasto-plastiques (Ottosen, Drucker-Prager), endommageables (Mazars, Laborderie), à fissuration diffuse (smeared crack model), de type gradient et de type énergétique (fictitious crack model (Hillerborg, 1984), crack band model). Dans certains cas, des modèles élasto-plastiques sont utilisés mais l’utilisateur de ce type de loi doit être très prudent une fois la phase élastique dépassée car il est possible d’obtenir des déformations qui ne reflètent pas du tout le comportement du béton.
  • Les modèles basés sur une approche discrète. En effet, la fissuration du béton engendre des discontinuités géométriques, et ces discontinuités sont intégrées entre les frontières des éléments connectés. Il existe les modèles discrets de Ngo, (Ngo, 1967), Bazant (Bazant, 1979), Blaauwendraad (Blaauwendraad, 1981), les modèles à éléments finis enrichis (E-FEM, X-FEM), (Belytschko, 1999) et également les modèles de type lattice.

Mais dans la littérature, les deux catégories de modèles les plus utilisées sont les modèles basés sur la mécanique de l'endommagement et les modèles considérant la présence de discontinuité explicitement.

Bien que la composante élastique ne soit pas exactement linéaire ([Baron] p.276 ; [Sargin 1971]), on effectue le plus souvent les calculs en adoptant l’élasticité linéaire isotrope (loi de Hooke) pour représenter la phase élastique du comportement du matériau béton. On adopte pour le module d’Young et le coefficient de Poisson des valeurs obtenues après caractérisation du béton. Dans le cas où il n’est pas possible de faire d’autres caractérisations mécaniques sur le matériau béton pour déterminer ses propriétés mécaniques actuelles, certaines de ses propriétés comme la résistance en compression peuvent être déduite à partir de lois-modèles que l’on trouve dans des codes-modèles comme l’Eurocode 2, par exemple.

Note - Les bureaux d’étude ne s’écartent de ce modèle que pour des études de détail, sur une partie de la structure, où des phénomènes mécaniques particuliers se sont manifestés : fissuration, réaction de gonflement, hétérogénéité locale. Il existe en effet des éléments finis de type « élément fissuré », des modèles d’endommagement, des lois viscoplastiques ou enfin des lois couplées poro-mécaniques ; mais il faut disposer, pour de telles études, de données issues de l’auscultation de l’ouvrage ou de prélèvements de béton pour s’orienter vers ces lois de comportement particulières, pour lesquelles il faudra de toute façon fixer plusieurs coefficients ou modules dans la préparation du calcul.


Modèle élastoplastique
 – De manière simplifiée, le béton est souvent modélisé comme un matériau élastoplastique homogène et isotrope, ce qui est, bien sûr acceptable tant que le béton ne fissure pas.

Vis-à-vis d’une loi élastoplastique, la fissuration du béton n’est pas directement modélisée : les zones de fissuration sont caractérisées par des déformations anélastiques importantes (>1-2‰) et un état de contrainte figé entre la résistance en traction ft et la résistance en compression fc du béton.

Comme on l’a vu au chapitre 3, une loi de comportement élastoplastique est composée d’une loi élastique et d’une loi d’écoulement plastique, associée à un « critère de plasticité ».

Le comportement post-fissuration peut néanmoins être représenté de façon approchée par une courbe σ-ε descendante au-delà de la résistance en traction. Celle-ci peut couvrir les effets de tension softening (représentant le travail nécessaire pour ouvrir la fissure) et les effets de tension stiffening (contribution du béton entre les fissures, contraint par l’adhérence à l’armature).

Si le comportement décroissant matériel conduit à un comportement décroissant global, il faut prêter attention à d’éventuels effets de localisation : la taille des éléments finis va limiter la taille des zones anélastiques et la solution va dépendre du maillage utilisé. Différentes techniques numériques permettent de résoudre ou limiter ce problème.

Endommagement - Une loi d’endommagement est une loi qui rend compte d’un des principaux effets macroscopiques de la fissuration du béton : la perte de raideur du matériau. L’idée fondamentale est de renoncer à suivre individuellement les fissures (leur apparition et leur propagation), et de considérer que le béton d’une structure périt plutôt par multiplication de fissures dans les zones dégradées. Ce type de loi permet donc de décrire la diminution de la rigidité du matériau sous l’effet de la création de micro-fissures dans le béton. Cette perte de rigidité est mesurée par une variable interne appelée endommagement, notée D, qui évolue de 0 (matériau sain) à 1 (matériau totalement endommagé) ; cette variable est généralement scalaire.

Afin de représenter au mieux le comportement du béton, les lois d’endommagement tiennent compte du comportement post-pic adoucissant. Cela permet de déterminer la contrainte de la manière suivante :

avec 

Le mérite de cette approche est de considérer le béton comme un « milieu continu », cadre auquel la méthode des éléments finis est parfaitement adaptée.

Dans le cadre du comportement du béton, les deux principales familles de modèles d’endommagement sont les modèles anisotropes et isotropes. L'isotropie caractérise l’invariance des propriétés physiques du béton quelle que soit la direction. Alors que l’anisotropie est dépendante de la direction. Une loi anisotrope présente des réponses différences dues au chargement selon son orientation.

Un des modèles d’endommagement les plus utilisés dans le monde industriel et de la recherche est le modèle de Mazars [Mazars, 1984]. C’est certainement le premier modèle d’endommagement du béton fonctionnant de manière correcte.

Les principales difficultés posées par les modèles d’endommagement sont :

  • une dépendance a priori des résultats au maillage : il convient en principe de démontrer que le mécanisme de ruine obtenu par ce type de modèle est indépendant de la finesse de maillage, au moins à partir d’un certain seuil. C’est d’ailleurs cette dépendance qui a amené au développement de méthodes de régularisation.
  • l’absence de solution analytique dans des cas même simples.

Plusieurs classes de régularisation existent dont la régularisation non locale et la régularisation par l’énergie de fissuration [Hillerborg, 1976] (qui ne résout que partiellement le problème). Parmi les méthodes non locales, il est possible de citer les méthodes intégrales [Pijaudier et al., 1987], [Giry et al., 2011] ou les méthodes à gradient (gradient de déformations ou gradients de variables internes [De Borst et al., 1992], [Peerlings et al., 1996], [Nedjar, 2005], [Lorentz2017]. Ces méthodes nécessitent d’utiliser des maillages relativement fins ce qui les rend généralement très coûteuses en temps de calculs.


Effets différés et redistribution des contraintes
 - Dans l’étude du comportement d’un ouvrage en béton sur les semaines suivant le coulage du béton, mais surtout sur le long terme, il est nécessaire de tenir compte des effets différés tels que le retrait et le fluage.

Ces phénomènes propres au béton peuvent, en principe, être modélisés en adoptant une loi de comportement visco-élastoplastique (loi de Bingham), ou parfois même seulement viscoélastique (fluage « scientifique », [Eymard]) : cette approche est généralement mise en œuvre par les laboratoires de recherche, pour analyser les essais sur le matériau. Mais, dans le cas d’une modélisation plus fine tenant de ces effets différés, il est nécessaire d’intégrer des phénomènes tels que le séchage et la fissuration en plus du retrait et du fluage. Car tous ces phénomènes interagissent entre eux, d’où leur importance dans une modélisation numérique.

La déformation engendrée par le retrait est une déformation induite par le séchage du béton provenant des effets de l’environnement. Ce retrait induit une déformation différentielle du béton autrement dit des contraintes plus élevées aux zones à fort séchage. Ce phénomène engendre des contraintes de traction en surface entraînant la fissuration et des contraintes de compression en cœur.

Concernant le fluage, sa déformation est généralement décomposée en deux déformations, l’une due au fluage propre et l’autre au fluage de dessiccation. Cela vient du fait que le fluage dépend fondamentalement de l’humidité relative.

Dans l’Eurocode 2, il est possible de déterminer la déformation du béton due aux effets différés (sans chargement extérieur). Pour ce faire, il faut calculer la déformation du au retrait endogène (conséquence de l’humidité interne) et la déformation de dessiccation (conséquence du séchage et de la taille de l’élément structurel). Et la déformation de fluage sous contrainte de compression σc est égale, selon le §3.1.9 de l’EN1992-1-1 à :

où Ec est le module tangent (égal à 1,05 Ecm) et φ est le coefficient de fluage.

L’Eurocode 2 (EN1992-1-1) donne une méthode simplifiée de calcul de φ(∞,t0). L'annexe B donne une méthode plus complète permettant d'estimer φ(t,t0) ainsi que l'évolution du retrait (voir également l'annexe B de l'EN1922-2). De plus, il est important de retenir que le calcul des déformations différées fait intervenir le type de ciment.


2) Acier de charpente

Le comportement de l’acier est beaucoup plus simple que celui du béton, et cela pour plusieurs raisons : c’est un matériau isotrope, présentant des résistances et des modules identiques en traction et en compression. De plus, il fait l’objet de contrôles industriels qui assurent son homogénéité.

Modèles élastiques - Bien que l’acier se comporte essentiellement comme un matériau élastoplastique isotrope (Eurocode 3, partie 1.1, §5.4.3), les modélisations usuelles des structures acier, ou composites comportant de l’acier, adoptent pour ce dernier un comportement élastique linéaire. L’on vérifie a posteriori que la contrainte dans l’acier est inférieure à la limite élastique fy. Dans le cas d’éléments de type poutre ou poteau, les normes autorisent, si les sections présentent des dimensions assurant une ductilité locale suffisante, à dépasser dans l’analyse la résistance élastique et à considérer des moments résistants basés sur une distribution plastique des contraintes.

Élastoplasticité et écrouissage - La théorie de l’élastoplasticité s’est en fait développée à partir de l’étude des alliages acier, notamment en vue de prédire les efforts de laminage et de forgeage ([Hill], [Nadai]). Pour l’acier, le modèle usuel se compose de la loi de Hooke pour la déformation élastique, et du critère de plasticité dit « critère de von Mises », comme y invite d’ailleurs l’Eurocode 3 (partie 2, § 7.3 ; partie 1.5 §10 et partie 1.7 §5.2.3.2) : pour ce critère, il suffit de donner la contrainte de limite élastique fy de l’acier.


Du point de vue du calcul par éléments finis, il se pose encore la question de la prise en compte de l’écrouissage du matériau, c’est-à-dire de son durcissement dans la phase purement plastique. Cet aspect est notamment abordé à l’annexe C6, partie 1.5 de l’Eurocode 3. Dans un souci de vérification du modèle, il est sans doute préférable d’effectuer un premier calcul avec un modèle sans écrouissage : en effet, on peut alors vérifier la qualité des contraintes obtenues dans les zones plastiques en affichant les isovaleurs de la contrainte de von Mises (les zones plastifiées doivent être à peu près monochromes). Il est cependant à noter que l’usage d’une loi élastique parfaitement plastique, sans écrouissage, peut entraîner des problèmes de convergence de l’analyse non linéaire. En effet, les zones plastifiées présentent un module tangent nul et n’ont donc pas de raideur. La prise en compte de l’écrouissage permet de stabiliser la résolution numérique.

En cas d’écrouissage jusqu’à une contrainte de valeur fy+X, les décharges et recharges éventuelles se font élastiquement, avec un module égal au module initial. fy+X devient donc la limite d’élasticité apparente du matériau. Par ailleurs, le comportement de l’acier est proche du modèle d’écrouissage cinématique. En première approche, l’on peut considérer qu’un écrouissage en traction jusqu’à une contrainte (fy+X) entraîne une diminution de la limite d’élasticité en compression à la valeur fy-X, et vice versa.

Le modèle le plus courant pour l’acier est cependant l’écrouissage isotrope, dépendant de la déformation plastique cumulée :


Dans ce cas, il n’y a plus de distinction entre traction et compression et un écrouissage dans un sens de sollicitation entraîne une augmentation de la limite d’élasticité en compression. Cette hypothèse n’est donc valable que si le chargement considéré est monotone, et non cyclique.

Dans presque tous les codes de calcul, on peut afficher les isovaleurs de la déformation plastique cumulée : ainsi, il est possible de vérifier, dans les zones plastiques, que la contrainte de von Mises se répartit conformément à cette déformation.

Dans le cas particulier de l’acier inoxydable, l’hypothèse d’un comportement linéaire jusqu’à la limite d’élasticité conventionnelle fy, correspondant à une déformation plastique de 0,2 %, n’est pas respectée, et, comme le précise l’Eurocode 3-1-4, §4.2., il convient de prendre en compte les effets du comportement non linéaire en contrainte-déformation dans le calcul des flèches. La loi matérielle est de type Ramberg-Osgood, et des indications pour sa modélisation sont données par l’annexe C de cette partie de l’Eurocode.

3) Acier de précontrainte

Du point de vue du comportement, les aciers de câbles (haubans, suspentes) et les aciers de précontrainte ne diffèrent des aciers de charpente que par une limite élastique environ trois fois plus élevée (de 1680 à 2140 MPa en France), justifiée d’ailleurs par des sollicitations en service généralement supérieures à 1000 MPa. On peut donc, en première intention, modéliser les armatures de précontrainte comme des éléments à comportement élastique linéaire, ou à comportement élasto-plastique de type von Mises.

Les aciers de précontrainte sont tous laminés à chaud, et il y a donc lieu, en phase plastique, de prendre en considération une loi d’écrouissage linéaire isotrope. Le coefficient d’écrouissage h peut être calculé à partir des caractéristiques garanties par l’armaturier :

où Rm est la contrainte à rupture, fp0,1 la limite conventionnelle d’élasticité, E le module d’Young et Agt l’allongement à rupture.


Une conséquence directe de ce niveau de contrainte élevé est l’amorçage de mécanismes de relaxation. La relaxation des aciers est un effet différé (dépendant essentiellement du temps écoulé depuis le chargement), donc non-élastique. Elle se traduit, pour une barre ou un câble soumis à un allongement constant, par une baisse progressive de la contrainte. Ce mécanisme ne s’amorce, à la température ambiante, que pour un allongement dépassant environ 60 % de la limite élastique (soit environ 1000 MPa). La relaxation augmente sensiblement avec la température.


Comme pour le fluage des bétons, on peut modéliser la relaxation avec une loi de comportement visco-élastique linéaire ; toutefois, cette approche est plutôt réservée aux travaux de recherche. En règle générale, il n’est pas courant de modéliser les aciers de précontrainte : on se borne souvent à introduire l’action de la précontrainte par une distribution de forces réparties à travers les éléments en béton. Cependant, pour des analyses de détail, où l’on s’intéresse à l’interaction entre les câbles et le coulis d’enrobage, on peut prendre en compte la relaxation par un calcul incrémental en temps, où les pertes de relaxation sont introduites comme des contraintes initiales.

4) Acier passif

Action mécanique normale des aciers passifs - Dans les calculs par éléments finis, les aciers passifs sont fréquemment modélisés par des éléments linéaires, c’est-à-dire unidimensionnels, de type « barre ». Suivant les principes des règlements de béton armé historiques, on considère en effet que les armatures agissent principalement en reprenant les tensions et compressions du béton suivant leur axe propre, et c’est précisément ce que permet le modèle d’élément fini de barre. L’accrochage aux nœuds des éléments massifs représentant le béton permet, au demeurant, d’éviter les mécanismes propres aux assemblages de barres seules.

Dans la plupart des codes de calcul, les éléments de barre sont considérés, par défaut, comme des éléments à comportement élastique linéaire. Étant donné le caractère unidimensionnel de ces éléments, on a donc une loi du type N = E A u/L, où u est le déplacement axial, et N l’effort normal aux nœuds. Il faut vérifier a posteriori que la contrainte axiale N/A reste comprise entre ±fy.

Adhérence - Ce qui précède ne vaut, bien sûr, qu’à la condition qu’il y ait adhérence continue entre les armatures passives et le béton.

5) Bois

Le bois est un matériau qui présente quelques singularités : il n'est pas homogène et ceci à différentes échelles, il n'est pas isotrope, il ne présente pas de comportement symétrique, il peut présenter des ruptures ductiles ou fragiles selon la sollicitation et son orientation. Le bois est sensible à l'humidité, ce qui a un impact sur ses caractéristiques dimensionnelles et sur sa raideur et sa résistance. La durée de chargement a un impact important sur la résistance et la déformée des éléments en bois ; les variations d'humidité peuvent accélérer ces déformées (mécanosorption). Sur la base de ce constat, on peut asseoir un raisonnement pour aborder une modélisation par la méthode des éléments finis.

Homogénéisation - L'hypothèse doit être validée sur le volume élémentaire représentatif, a minima, le volume du plus petit élément fini. En sachant que, selon les essences, les cernes de croissance peuvent dépasser 1 cm, il deviendra problématique d'assumer cette hypothèse d'homogénéisation à proximité des organes d'assemblages (pointes, broches, boulons…) qui présentent un diamètre de cet ordre de grandeur, voire plus faible.

La présence de nœuds est rarement prise en compte dans les modélisations de structure ou d'éléments de structure.

Orthotropie - Le bois présente une structure et des caractéristiques fonction de trois directions, la direction longitudinale – l'axe du tronc – le fil du bois, les directions radiales et tangentielles, perpendiculaires à la direction longitudinale. Ces dernières s'inscrivent dans un plan, souvent celui des sections droites de poutres, sur lequel apparaît plus ou moins nettement les cernes de croissance. Le repère d'orthotropie est donc un repère quasi cylindrique, alors que le repère des éléments s'inscrit plutôt dans un référentiel cartésien. La représentation de cette orthotropie, quand elle est prise en compte, se limite dans la grande majorité des cas à une hypothèse d'isotropie transverse (axes radial et tangentiel à caractéristiques identiques) dans un référentiel cartésien. La pente de fil associée à la présence de nœuds n'est en général pas prise en compte dans le calcul des éléments. Par contre, elle est modélisée dans le calcul des assemblages, notamment pour les assemblages par contact.


La matrice de souplesse (inverse de la matrice de rigidité) peut être définie comme suit pour un cas plan orthotrope :

Remarque : L’hypothèse d’orthotropie et de symétrie de Sij réduit le nombre de termes indépendants de 36 (cas 3D le plus général) à 9 termes.

Bois – hypothèse orthotrope


Modélisation élastique
 - Il suffit de disposer de la matrice de comportement isotrope transverse pour une modélisation 2D ou 3D. Les modules entre les directions longitudinale et radiale ou tangentielles peuvent présenter des rapports de l'ordre de 20. Ils évoluent en fonctions de la durée de chargement (fluage), de l'humidité du bois à la mise en œuvre et de son évolution dans le temps (environnement). La représentativité du résultat MEF dépendra de la pertinence des paramètres pris en compte.

Modélisation plastique, critères de rupture - Pour une matériau isotrope transverse non symétrique, on peut s'orienter vers des critères de Hill, de Tsai… en sachant qu'il sera nécessaire de décrire la rupture fragile, en traction perpendiculaire et en cisaillement. La grande variabilité des résistances rend le calage des paramètres des critères délicat. La succession de cernes de croissance, ou celle des lamelles assemblées, avec des caractéristiques mécaniques différentes, peut rendre un modèle homogène difficile à caler en termes de résistance. En effet, la résistance et la rigidité sont fortement corrélées pour le bois et des effets « système » apparaissent rapidement en termes de résistance d'élément. Ainsi, on dispose de limites de résistance en traction longitudinale, en compression longitudinale, en flexion…


Enfin, les structures en bois sont particulièrement sensibles au comportement de leurs assemblages. Ceci peut présenter l'avantage de ne pas avoir à modéliser finement les éléments, mais à s'intéresser plus précisément aux zones d'assemblages. Par contre, apparaissent les problématiques de contacts, de matériaux multiples, de la plastification de certains, de rupture fragile d'autres, des limites d'homogénéisation rappelées précédemment.

Il apparaît clairement que l'effort de modélisation est lié à l'échelle d'investigation, ou au stade du projet. Sous réserve de prise en compte des singularités listées ici, la modélisation du bois au sein d'une structure peut être menée comme pour un autre matériau.

Modélisation des effets différés et de l’interaction avec les phénomènes hydriques - Le bois est un matériau hygroscopique sensible au changement d’humidité relative de l’air. Il est par ailleurs sujet à des effets de fluage sous contrainte. Si l’on désire modéliser ces phénomènes, l’on peut adopter un modèle de type viscoélastique vieillissant en environnement variable, en accord avec les principes de la thermodynamique. Le modèle de comportement viscoélastique de type Kelvin-Voigt généralisé, caractérisé par des paramètres rhéologiques vieillissants et dépendant du niveau et de l’histoire de l’humidité peut être associé au modèle de comportement non-viscoélastique de Ranta-Maunus pour caractériser le retrait-gonflement et la mécano-sorption.

$translationBooks

D.2 Les catégories d'éléments d'ouvrage

D.2 Les catégories d'éléments d'ouvrage

On présente ici les particularités de ces éléments eu égard à un calcul aux éléments finis.

1) Éléments de Béton Armé

Prise en compte des phénomènes différés - Le plus souvent, l’ingénieur de projet s’intéresse aux effets de redistribution des contraintes et aux déformations différées qui accompagnent le vieillissement des matériaux.

Pour y parvenir, le calcul statique par incréments de temps est particulièrement bien adapté et largement suffisant. Le retrait et le fluage du béton, par exemple, ne dépendent que du temps écoulé depuis le remplissage du coffrage. Dans un calcul incrémental, ils seront introduits comme des déformations imposées en chaque nœud du maillage. Il est possible de calculer a priori, pour chaque pas de temps ou incrément, la carte des déformations de retrait et de fluage : le solveur éléments finis intègre cet état initial dans la recherche de l’équilibre d’un matériau élastoplastique.

Il convient néanmoins de prêter attention à l’interaction entre les effets différés et les phases de construction, comme cela sera expliqué par la suite.

Prise en compte de la fissuration - Lorsque le béton se fissure en traction, il se développe des fissures débouchant vers le parement le plus proche, ainsi que des fissures le long de l’interface acier béton ([Goto], figure ci-dessous).

Fissuration interne du béton armé

Lorsque l’on considère une bielle tendue de béton armé fissuré comme un milieu continu homogénéisé, la loi de comportement N(ε) (où N est l’effort normal de traction) suit l’allure de la figure ci-dessous.

Modèle schématique de l'effet de tension stiffening

Le Model Code 1990 de la CEB-FIP en donne une expression analytique. Cette relation peut être utilisée en modélisant la section de béton fissurée comme un élément de barre dont on actualiserait le module de rigidité par incréments (suivant la valeur de l’effort normal ou de l’allongement, puisque la relation est inversible). On peut aussi l’utiliser comme loi de comportement d’une fibre d’un élément multifibre.

Béton confiné - Pour les zones de béton confinées (par des cadres d’armature passives, par exemple), il est possible de prendre en compte la résistance résiduelle post-écrasement, comme l’indique l’EN1992-1, en utilisant la loi de Sargin :

 avec εc1 la valeur au pic telle que 

 avec 

Pour 

L’EN1992-1 au §3.1.9 propose également un incrément de résistance et de déformation du béton quand il est soumis à un contrainte de confinement de valeur σ23 tel que :

Cette loi n’est toutefois qu’unidimensionnelle, et ne vaut que sous chargement monotone. On ne peut donc s’en servir avec la méthode des éléments finis qu’en modélisant le béton confiné par un élément de barre (donc on modifierait incrémentalement le module), ou comme loi de comportement d’une fibre d’un élément de poutre multi fibres.


2) Éléments en Béton Précontraint en Pré-tension

La précontrainte par pré-tension, caractéristique notamment des produits en béton préfabriqués industriellement (poutres, prédalles, dalles alvéolées, poutrelles…), consiste à venir tendre des câbles (fils, torons ou barres) dans des bancs de fabrication puis de venir couler le béton avant de détendre les câbles lorsqu’une résistance minimale du béton est atteinte (appelée résistance du béton au relâchement). L’intensité de la tension des câbles (qui ne doit pas dépasser la force de précontrainte maximale autorisée par les codes de calculs), le nombre de câbles et la résistance du béton sont ajustés en fonction des charges que doit reprendre le plancher ou l’élément précontraint.


Lors de cette détention, des pertes instantanées, de l’ordre de 8% pour une prédalle précontrainte par exemple, doivent être prises en compte (pertes dues à la rentrée d’ancrage, à la relaxation des armatures de précontrainte pendant la période entre la mise en tension des armatures et le transfert de la précontrainte, au raccourcissement élastique du béton sous l’effort de compression imposé par la précontrainte) lors d’un dimensionnement en phase provisoire de chantier avec une longueur de transmission de la précontrainte à considérer en partant de l’about de l’élément préfabriqué en béton.

A plus long terme, des pertes différées dues au retrait du béton, à la relaxation de l’acier ou au fluage du béton pour atteindre au final une perte totale de 20 % par exemple pour une prédalle précontrainte, doivent également être considérées.


Dans de nombreux codes de calculs aux éléments finis, la précontrainte peut être intégrée dans des éléments finis de type poutre représentant les câbles liaisonnés à des éléments finis volumiques représentant le béton. Suivant la zone d’étude, il peut être nécessaire par exemple dans les zones d’abouts de considérer la longueur de transmission de la précontrainte dans les câbles. La distribution de la force réelle de précontrainte (avec prise en compte des pertes instantanées ou différées en fonction du moment dans la durée de vie du produit où l’on souhaite faire le calcul/vérification) est alors variable le long de cette longueur de transfert. Une distribution linéaire de cette force de précontrainte est autorisée dans la plupart des codes de calculs et reste dans la plupart des cas sécuritaire pour le dimensionnement par rapport à la distribution parabolique plus réaliste.


Pour des raisons de complexité et de besoin (dimensionnement limité à des déformations élastiques), la modélisation de ces éléments est réalisée la plupart du temps avec des hypothèses linéaires (lois de comportement du béton et de l’acier élastique linéaire, contact parfait entre le béton et l’acier…). Pour des études fines, des hypothèses non linéaires peuvent toutefois être mises en place selon le besoin comme par exemple des lois de comportement de type endommagement pour le béton, de type élasto-plastique pour l’acier et l’introduction d’éléments d’interface acier-béton.


3) Éléments en Béton Précontraint en Post-tension

Pertes élastiques - La mise en tension des câbles de précontrainte par post-tension s ’accompagne de pertes instantanées : frottement, recul d’ancrage et perte par allongement élastique.

Effets différés : fluage, retrait et relaxation - La prise en compte des effets différés s’effectue par un calcul incrémental. Le fluage et le retrait peuvent être introduits comme des déformations volumiques incrémentales, données en chaque nœud du maillage.


4) Éléments de charpente métallique

Choix du type d’analyse

De par leur fort élancement, les charpentes métalliques sont très déformables. Il en résulte que l’hypothèse traditionnelle de réaliser l’équilibre des efforts en configuration initiale n’est pas toujours valable ; il faut alors établir les distributions d’efforts internes en configuration déformée.

La sensibilité à ces effets non linéaires, dits parfois du second ordre, est jugée au travers du multiplicateur critique αcr, le multiplicateur des charges conduisant à l’instabilité eulérienne de la structure. Dans la version actuelle de l’Eurocode 3 :

  • si αcr > 10, les effets non linéaires peuvent être négligés. Si l’analyse structurelle globale intègre la plasticité des éléments, la valeur limite de αcr est augmentée jusqu’à 15 ;
  • si 4 < αcr < 10, ils doivent être pris en compte. Ceci peut cependant se faire par une analyse élastique linéaire classique, en amplifiant les efforts transversaux ;
  • si αcr < 4, l’analyse non linéaire est obligatoire.

Dans les deux derniers cas, les imperfections géométriques globales doivent être considérées, ainsi que les imperfections d’éléments si elles influent sur le comportement global.


Les imperfections d’éléments regroupent :

  • les imperfections géométriques : imperfections transversales et en torsion ;
  • les imperfections matérielles : les éléments laminés, ou reconstitués par soudage, présentent des distributions de contraintes résiduelles autoéquilibrées créées par leur processus de fabrication.

Ces dernières peuvent être représentées par une imperfection géométrique équivalente, dont la valeur peut être trouvée dans les normes en vigueur.

Ces imperfections doivent être intégrées pour toute analyse d’un élément de charpente intégrant les effets non linéaires.

Les modèles de type structure ou RdM (barres, poutres, plaques et coques) sont en principe idéalement adaptés aux calculs de charpente métallique.


Structures filaires

Analyse élastique

En vue de la modélisation, il y a lieu toutefois d’analyser précisément :

  • la nature des liaisons entre les différents éléments de charpente,
  • et le mode de transfert des charges de chaque pièce aux autres.

Analyse plastique

Lorsque la ductilité sectionnelle et/ou la ductilité des assemblages est assurée, il est possible de réaliser des analyses structurelles intégrant la plasticité. Différentes méthodes peuvent être utilisées (EC3, §5.4.3) :

  • analyse élastique-plastique, où les zones plastifiées sont modélisées comme des rotules plastiques ;
  • analyse plastique non linéaire, prenant en compte la plastification partielle des barres le long des zones plastiques ;
  • analyse rigide-plastique, dans laquelle le comportement élastique des barres entre les rotules est négligé.


Torsion

Les éléments en acier sont rarement massifs ; ils sont composés de parois minces pour constituer des profils, ouverts ou fermés. En particulier dans le premier cas, la réponse aux sollicitations en torsion se fait à la fois en torsion uniforme, dite de Saint Venant, et en torsion non uniforme, entraînant un gauchissement de la section. Ce dernier phénomène n’est habituellement pas pris en compte dans les logiciels commerciaux, alors qu’il peut avoir un effet déterminant dans la réponse des structures. Dans ce cas, deux solutions sont envisageables :

  • procéder à une modélisation surfacique de l’élément : cette solution est inapplicable si l’analyse porte sur une structure comportant plus de quelques éléments ;
  • si l’on est dans une situation similaire à celle de la poutre en I, où la torsion non uniforme peut se représenter par la flexion alternée des semelles (appelée biflexion, Figure), procéder à une modélisation bifilaire de l’élément : les deux semelles sont représentées par deux éléments distincts, reliés par des éléments transversaux représentant l’âme. De ce fait, la flexion spécifique de chaque semelle, et partant la torsion non uniforme, est représentée. Ce cas se présente notamment pour les ponts mixtes bipoutres. Une application est présentée au chapitre 3.

Décomposition de la torsion en torsion uniforme et non uniforme : hypothèse simplifiée de biflexion


Éléments bi ou tridimensionnels

L’analyse élastique linéaire des éléments aciers bi ou tridimensionnels ne présente pas de problèmes spécifiques, et les règles générales s’appliquent.

Par contre l’analyse non linéaire peut être requise pour étudier les phénomènes d’instabilité. En effet, les structures composées de plaques d’acier, qu’elles soient planes ou courbes, sont sujettes à des phénomènes de voilement.

Dans le cas des éléments plans, appelés plaques, le voilement est un phénomène relativement stable : l’initiation du voilement de la plaque n’entraine pas la ruine, la charge maximale est atteinte après voilement. On parle de comportement post-critique.

Dans le cas des éléments courbes (coques), l’instabilité entraine la ruine immédiate de la structure, souvent brutale. D’un point de vue numérique, dans une analyse statique non linéaire, cela se traduit par une rapide décroissance de la charge au-delà du maximum.

Dans les deux cas, la charge maximale est fortement dépendante de la déformée initiale appliquée, tant en amplitude qu’en forme. L’amplitude est fixée par les normes. La forme est choisie habituellement affine au premier mode d’instabilité. Cependant ce choix n’est pas nécessairement le plus pénalisant. Il est conseillé de le compléter par des modes locaux lorsque la structure présente des panneaux de dimensions fortement différentes.

Par exemple, dans le cas d’un platelage orthotrope, la déformée affine au mode de voilement global doit être complétée de déformées affines au voilement des sous panneaux.

Dans le cas des structures de type coque, le problème est encore plus critique. Il est conseillé, une fois un premier calcul effectué selon les hypothèses ci-dessus, d’adopter dans un deuxième calcul une forme de déformée initiale affine à la déformée obtenue à la ruine.

5) Structures mixtes : acier-béton

L’alternative d’une construction mixte acier-béton est parfois préférée pour certains types de bâtiments industriels et pour des ponts de petites à moyennes portées (travée centrale < 100 m). L’alliance de ces deux matériaux en les faisant « travailler » dans leurs domaines de résistance (béton en compression et acier en traction) permet d’obtenir une structure légère et résistante. Afin d’aboutir à ce résultat, la connexion entre ces deux matériaux se doit d’être correctement dimensionnée. On distingue :

  • les dalles mixtes : dalle pleine + bac collaborant,
  • les poutres mixtes : dalle pleine ou mixte + profilé métallique + connecteurs,
  • les poteaux mixtes : profilé métallique + (remplissage béton ou enrobage béton).

Les analyses structurelles globales sont usuellement menées par analyse élastique en homogénéisant la section, ou en représentant de façon séparée les deux matériaux. Cette deuxième façon de procéder peut conduire à des difficultés dans le traitement des résultats, puisqu’il faut recalculer les efforts sur la section globalisée pour pouvoir appliquer les règlements en vigueur.


Lorsque l’on désire effectuer des analyses plus fines, étant donné la diversité des matériaux mis en collaboration dans ce type d’alliance, autant sur le plan géométrique que sur le plan du comportement en général non-linéaire, des modèles aux éléments finis 3D s’avèrent nécessaires pour mener des études locales incluant un traitement des diverses interfaces (par des éléments finis de contact par exemple). Pour des études à l’échelle de la structure ou des éléments, des modèles 2D relativement performants ont été développés cette dernière décennie notamment, ceux basés sur un découpage de la section transversale en fibres (modèle à fibres) pour permettre d’estimer par intégration numérique la rigidité en section.

Le raboutage des poutres de ponts mixtes constitue aussi un détail relativement complexe. Quelle que soit la solution de raboutage retenue (par couvre-joints, par chevêtre ou par diaphragme), les modèles 3D utilisant des éléments finis de type massifs sont préférés aux modèles 2D simplifiés mais requièrent des temps de calculs relativement importants.


Fissuration
 - Les poutres mixtes sont usuellement constituées d’un profil métallique solidarisé à un plancher ou un tablier en béton armé au moyen de connecteurs. Il en résulte que, dans les zones de moment positif, où la dalle est comprimée, la résistance et la raideur sont particulièrement importants, alors que dans les zones de moment négatif, la fissuration entraîne des caractéristiques mécaniques nettement moindres. Ceci ne peut être négligé dans l’analyse structurelle globale de la structure. Différents niveaux de modélisation sont admis dans les codes :

  • Approches forfaitaires : par exemple, l’Eurocode 4 préconise de prendre une longueur fissurée égale à 15 % de la portée, de part et d’autre des appuis. Il propose aussi une raideur forfaitaire pour les poteaux mixtes. L'eurocode 8 quant à lui préconise d’adopter une raideur moyenne sur toute la longueur de la poutre. Ces différences d’approche se justifient par la forme différente des diagrammes de moments sous sollicitation classique majoritairement gravitaire, et sous sollicitation sismique ;
  • Approches définissant une zone de fissuration par analyse des sollicitations enveloppes : l’Eurocode 4 préconise de considérer comme fissurée toute section où la contrainte dépasse de deux fois la résistance en traction moyenne sous l’enveloppe des sollicitations caractéristiques calculées en supposant la structure non fissurée, en adoptant un module de béton à long terme
  • Il est aussi loisible de recourir à des analyses non linéaires.


Connexion

Sauf en cas d’analyse non linéaire, il n’est pas utile de modéliser la connexion. Les normes actuelles permettent de prendre en compte l’effet d’une connexion partielle sur la résistance des éléments.

Dans le cadre de la modélisation, des éléments finis de type poutre (2D ou 3D) sont en général utilisés pour modéliser la connexion ponctuelle. Pour une hypothèse de connexion répartie, il existe des modèles appropriés dans la littérature.


Largeurs collaborantes

Les poutres acier sont connectées à des éléments particulièrement larges et le trainage de cisaillement peut entraîner une répartition non uniforme des contraintes sur la largeur de la dalle.

Dans les modélisations filaires, ce phénomène est habituellement pris en compte en adoptant une largeur réduite de dalle dans la modélisation, à contrainte constante puisque la modélisation poutre l’impose.

En toute rigueur, comme le trainage de cisaillement est lié à la transmission de l’effort rasant par la connexion à la dalle, et est donc dépendant de la forme du diagramme de moment, la largeur collaborante devrait varier de combinaison à combinaison. Cependant les normes autorisent d’adopter une largeur unique pour l’ensemble des calculs.

Cette variabilité des largeurs collaborantes sera prise en compte si une modélisation de la dalle est réalisée par des éléments de type surfacique.


Effets différés et retrait

Les effets différés du béton influent sur la répartition des raideurs dans la structure et donc sur la répartition des efforts et doivent être pris en compte. Il est à noter que l’Eurocode module la valeur du coefficient d‘équivalence acier-béton en fonction du type de chargement. La sollicitation de retrait induit une distribution de contrainte autoéquilibrée sur la section ; celle-ci aussi doit être prise en compte.


Analyses non linéaires

La modélisation non linéaire des structures mixtes adopte les hypothèses matérielles et d’imperfections géométriques utilisées pour les matériaux béton et acier seul.

Comme déjà dit, la connexion doit être modélisée avec sa propre raideur. Il est à noter que son traitement numérique reste délicat : habituellement, il est supposé que les parties béton et acier ne peuvent se désolidariser transversalement, et que seul un glissement longitudinal est possible. La formulation éléments finis d’une telle liaison peut entraîner des phénomènes de « locking » : le glissement est bloqué par un jeu de projection de la raideur transversale empêchant le soulèvement lorsque l’équilibre est établi en configuration déformée. Il convient donc, lorsque ce type d’élément est utilisé, de vérifier la cohérence des efforts apparaissant dans la connexion avec les efforts de traction/compression des éléments acier.


6) Haubans, suspentes et câbles porteurs


Introduction dans les calculs
 - Certains éléments de câbles peuvent être modélisés comme des barres de section équivalente : tel est le cas des suspentes verticales, ou des câbles de précontraintes, guidés dans leur gaine : la courbure de ces câbles à l'équilibre est pratiquement indépendante de leur masse linéique. Dans les autres cas, le poids propre des câbles tend à la courber : la direction verticale est une direction particulière qui intervient dans la raideur de ces éléments.

Lorsqu’à l’échelle de la portée du câble, sa courbure devient importante, celle-ci conditionne fortement les contraintes et les efforts qu’il peut transmettre au reste de l’ouvrage. Pour les haubans, par exemple, la raideur du câble dépend principalement de sa flèche, et de son allongement. La flèche dégrade la raideur du câble, et donc la relation tension/flèche du câble est par essence non-linéaire.

Les éléments de câble présents dans la plupart des codes s’inspirent du modèle de [Gimsing], qui suppose que le câble se déforme selon une courbe parabolique, hypothèse qui est valide dès lors que le rapport flèche/portée est inférieur à 1/12 environ. Ce modèle découple l'allongement élastique de la flexion du câble.


Effet de chaînette
 - Lorsque la flèche d'un câble est supérieure à 1/12 de sa portée, il n'est plus possible de découpler l'allongement élastique de la flexion, car le bras de levier de la tension devient prépondérant dans l'expression du moment fléchissant. Plusieurs codes de calcul proposent des éléments de caténaire pour rendre compte de cet effet géométrique. Ces éléments connectent deux nœuds du maillage, et un seul élément suffit pour modéliser le câble.

Déplacements nodaux d'un élément de chaînette

Pour ce type d'élément, les inconnues nodales sont le déplacement vertical (dans la direction de la gravité) et le déplacement horizontal (cf. figure 4) ; les réactions associées sont la variation de la composante horizontale de la tension du câble, et la réaction verticale aux ancrages. La relation entre les déplacements nodaux et ces réactions dépendant de la tension du câble, la recherche de l'équilibre est un problème non-linéaire, quoique la structure soit globalement élastique.

$translationBooks

D.3 Les phases de construction

D.3 Les phases de construction

Les études des phases de construction ont deux objectifs :

  • s’assurer de la stabilité de la structure dans les différents états transitoires conduisant à l’état définitif ;
  • calculer les effets du montage sur la distribution des efforts et sur la déformée de la structure en service.

Les effets du montage sont multiples. Ils sont liés :

  • à l’évolution du schéma statique en cours de construction. Par exemple, une travée de pont construite par deux encorbellements à partir de ses deux extrémités et finalement clavée en son centre présentera juste après le clavage un diagramme de moment sous charges permanentes qui s’annule en milieu de travée, bien différent de celui qui aurait été obtenu sans prise en compte du mode de montage (Figure 5) ;
  • à l’interaction des effets différés avec l’évolution du schéma statique. Dans l’exemple cité précédemment, après clavage, le moment en milieu de travée va croître sous l’effet du fluage ;
  • à l’évolution des sections au cours du temps. Par exemple, dans le cas des ouvrages mixtes acier-béton, le poids de la dalle est supporté par la structure métallique seule si la structure ne repose pas sur des cintres lors du coulage du béton;
  • à des réglages volontaires de la structure : dénivellations d’appui, réglages des haubans, vérinages en clef d’arc, …


La prise en compte de ces effets peut s’avérer relativement complexe et, dans les cas les plus difficiles, il peut être indispensable de recourir à un logiciel capable de modéliser pas à pas l’évolution de la structure.

Il est cependant souvent possible de procéder par une superposition de différentes analyses linéaires.

Effet des phases de construction : principe dans le cas du clavage d'une travée de pont construite en console

Les principales difficultés sont liées aux effets des déformations différées du béton. En effet, comment évaluer le moment à terme en milieu de travée, dans l’exemple présenté plus haut ? Dans le cas d’un clavage unique, il est possible d’utiliser la méthode dite « des coefficients » (Figure 6). Cette approche repose sur les points suivants :


Etat à terme = (E(t0, t1)/E(t0, t)) x Etat à terme non clavé + (1-(E(t0, t1)/E(t0, t))) x Etat à terme sans phasage

t0 étant le temps d’application de la charge, t1 le temps du clavage, t le temps considéré pour l’état à terme, et E(t0, t1) le module de béton pour l’obtention de la déformation du béton au temps t1 pour une contrainte appliquée en t0.

Cette méthode, dans le cas d’un clavage unique, restitue l’état à terme théorique exact. Elle est cependant difficilement extensible au cas de clavages multiples, et peut conduire alors à des résultats aberrants.

Cas unitaires utilisés pour la méthode des coefficients

Il est préférable d’extérioriser les effets des modifications du schéma statique de la façon suivante (Figure 7) :

  • cas 1 : calcul de l’état à terme, si le clavage n’était pas réalisé ;
  • cas 2 : calcul juste avant clavage, avec le module de béton adéquat ;
  • cas 3 : calcul de l’effet du clavage : appliquer à la structure un déplacement imposé au niveau du clavage, ramenant la valeur de la discontinuité (dans le cas de l’exemple, en rotation) à la valeur figée par le clavage ;
  • l’état à terme est la somme des cas 1 et 3.

Cette technique est plus facilement extensible aux cas où les modifications de la configuration statique des ouvrages sont nombreuses.

Cas unitaires utilisés pour la méthode de superposition

Il n’est pas inutile de rappeler que, sous les effets de déformation différée, le béton réagit avec un module de déformation apparent, appelé module de relaxation, plus faible que le module de fluage correspondant. Si le rapport classique entre le module de l’acier et du béton, intégrant les effets de fluage, est de l’ordre de 18, dans le cas d’un déplacement imposé, ce rapport passe à une valeur de 24. Ceci a tendance à rendre moins efficaces les réglages par dénivellation d’appui et par vérinage, lorsque ce dernier conduit à imposer une déformation à la structure.

Du point de vue des simulations éléments finis, le nombre important d’états intermédiaires à traiter multiplie le risque d’erreurs. Les vérifications doivent porter :

  • sur le respect des conditions limites dans les phases intermédiaires ;
  • sur le respect des déplacements gelés par les phases de construction dans la structure (par exemple, discontinuité de pente au niveau de clavages).

En études d’exécution, il faut aussi rappeler que le fluage réel du béton peut s’écarter fortement des formulations théoriques ; il faut donc construire le modèle de telle façon qu’il soit aisé de l’adapter pour restituer les déformations apparaissant dans les premières phases, et ainsi améliorer la prédiction des phases suivantes.

$translationBooks

Chapitre E - Post-traitements typiques du génie civil

Chapitre E - Post-traitements typiques du génie civil


Introduction

La variété des quantités d’intérêt analysées pour des ouvrages de génie civil est importante car elle découle de la grande diversité de ce qui peut être étudié compte tenu de :

  • la nature de l’ouvrage considéré et donc des exigences fonctionnelles qui lui sont associées (par exemple l’étanchéité pour un barrage, un réservoir, une enceinte de confinement de centrale nucléaire (ouverture de fissure, état de contrainte, déformations résiduelles, …),
  • les états limites considérés (ultimes ou de services, …),
  • la nature des cas de charge (dynamique, statique, différé, …),
  • les éléments structurels constitutifs (béton massif, armé, précontraint, construction métallique, bois, ouvrages en maçonnerie, …).

Par ailleurs, les quantités d’intérêt requises peuvent être directement ou non accessibles à l’issue d’un calcul EF. Elles peuvent être construites à partir de produits du calcul qui pourront permettre par post-traitement de fournir la quantité d’intérêt qui sera ensuite comparée à un critère ou qui sera ensuite une donnée pour la suite du projet (sections d’armatures, orientations de barres de ferraillages, …).

Les Eurocodes permettent de savoir quelles sont les grandeurs à analyser compte tenu des différents cas cités plus haut. En revanche, il n’est pas établi dans les codes et normes comment avoir accès à ces grandeurs. Ce chapitre permet de fournir les clés permettant de savoir comment y accéder et les pièges à éviter.

E.1 Généralités

E.1 Généralités

E.2 Grandeurs en dynamique

E.2 Grandeurs en dynamique

E.3 Cas spécifique du béton armé

E.3 Cas spécifique du béton armé

$translationBooks

E.1 Généralités

E.1 Généralités

Les différentes grandeurs post-traitées peuvent différer très significativement en fonction de la finalité de l’analyse faite par éléments finis.

Typiquement, les grandeurs étudiées au titre d’études de prédimensionnement (esquisse ou Avant Projet Sommaire) sont très différentes de ce que l’on va chercher à obtenir en Avant Projet Détaillé et a fortiori en études d’exécution. Lorsqu’il s’agit d’analyser le comportement d’un ouvrage existant pour quantifier un éventuel risque structurel les quantités d’intérêt sont encore différentes.

On serait tenté de considérer :

  • Phases de prédimensionnement (esquisse, APS, …) : grandeurs de type déplacements, flèches rapportées à des portées, etc …,
  • Phases de dimensionnement : densités de ferraillage pour des plaques, sections de ferraillage pour des poutres et poteaux, ouverture de fissure, …
  • Analyse structurelle : états de contraintes (contraintes principales dans le béton: directions, signes, tri-axialité et valeurs, comparaison à des critères, par exemple Ottosen, Rankine, Drucker-Prager si l’analyse est élastique), déformations principales directions, valeurs et signes, cartes d’endommagement, contraintes dans les aciers, déformations plastiques, …

Enfin, il importe de garder à l’esprit que toutes les quantités d’intérêt ne sont pas accessibles pour tous les différents types d’éléments et modèle de comportement. En effet, le type d’élément conditionne la nature des degrés de liberté : on n’a pas accès à des rotations aux nœuds et pas de moments vec des éléments volumiques de base. De la même façon, un modèle de plasticité ne peut fournir de valeurs d’endommagement. Ces éléments peuvent paraître triviaux mais ils sont suffisamment souvent oubliés par les utilisateurs des codes de calcul qu’il nous paraissait utile de le rappeler ici.


1. Contraintes et déformations

Comme cela est exposé dans le chapitre « Généralités », les champs de contraintes et déformations ne sont pas des produits de calcul directs d’une résolution d’un problème de mécanique par un solveur aux éléments finis. Ils sont déduits des déplacements par dérivations des champs interpolés de déplacements. On note en conservant les notations du chapitre généralités :

En élasticité linéaire, la relation contrainte – déformation s'écrit (avec H la matrice de Hooke) :

Dans ce genre de démarche, les champs de contrainte et de déformation évoluent selon des fonctions d’interpolation de degré inférieur à celle des déplacements.

En analyse non linéaire les contraintes et déformations ne peuvent pas être ainsi évaluées par un produit de la matrice de Hooke et de la dérivée des déplacements nodaux et s’appuient sur les valeurs aux points de Gauss pour un élément en faisant appel aux fonctions de forme pour les extrapoler sur le reste de l’élément et aux nœuds. Cette démarche permet en outre de faire en sorte que les contraintes évoluent selon le même degré que les fonctions d’interpolation et non leurs dérivées.

A titre d’exemple, sur un élément qui a n nœuds, npg points de Gauss situés aux coordonnées ξnpg, munis de fonctions de forme Ni, on retient la minimisation au sens des moindres carrés entre le champ interpolé évalué à partir des valeurs nodales recherchées et les valeurs gaussiennes connues.


Construction des valeurs nodales élémentaires à partir des valeurs aux points de Gauss en 1D

Soit donc il s’agit de minimiser la fonctionnelle :

Soit pour chaque nœud i parmi les n nœuds de l’élément :

Qui peut se mettre sous une forme matricielle dont les matrices sont connues une fois pour toute pour les éléments isoparamétriques de référence :

Et donc, on obtient directement les valeurs nodales de contrainte :


Il peut donc exister plusieurs valeurs nodales de déformation ou de contrainte pour un nœud commun à plusieurs éléments. Les valeurs nodales doivent être déduites à partir de ces valeurs.

Discontinuité des valeurs nodales élémentaires en 1D


Il existe diverses méthodes permettant de reconstituer la continuité des déformations ou des contraintes (si les matériaux sont les mêmes entre deux éléments contigus) en ayant une seule valeur nodale de déformation (ou de contrainte). Il s’agit de lissage des champs de contrainte et de déformation par certains codes de calculs sur la globalité de la géométrie (par minimisation au sens des moindres carrés) ou par des moyennes des valeurs provenant pour un nœud des éléments auxquels il appartient tel que cela figure dans la figure 3 pour un cas 1D.

Exemple de construction d’une valeur nodale moyenne des valeurs nodales élémentaires en 1D

Pour des raisons de visualisations on peut préférer de faire apparaître des grandeurs élémentaires lissées. Cependant, lorsqu’il s’agit de rechercher de la précision les grandeurs calculées et déduites au point de Gauss seront à privilégier.


2. Passage des contraintes de la mécanique des milieux continus aux efforts de la mécanique des structures

Les méthodes automatisées de calcul de ferraillage s’appuient sur des algorithmes en flexion composée sur base de torseurs de type plaque et coque à 8 composantes ou en flexion composée déviée sur base torseurs de type poutre à 6 composantes.

Certains ouvrages peuvent nécessiter pour diverses raisons de prendre en compte une modélisation volumique de l’ouvrage (barrage, structures de types réservoirs précontraints : réservoir de Gaz Naturel Liquéfié, enceintes de confinement). Pour décliner une démarche calculatoire sur base de torseurs à 8 composantes, il est nécessaire de reconstruire des efforts internes de mécanique des structures à partir des champs de contraintes le long de segments comme cela est exposé figure ci-dessous.

Segment le long duquel les contraintes servent de référence pour la reconstruction d’un torseur coque

On considère le tenseur des contraintes exprimé le long du segment de la figure précédente :

En toute rigueur, pour une plaque mince satisfaisant les hypothèses de Kirschhof-Loeve, on considère le torseur coque à 10 composantes que l’on associe aux contraintes suivant le principe d’équivalence :

Cette intégration continue est en réalité discrète le long du segment et typiquement peut se réaliser à parti des valeurs des contraintes exprimées au point de Gauss ou aux nœuds.

Diverses techniques d’intégration peuvent être utilisées (Trapèze, Gauss, Lobatto, etc …).

On conservera par la suite et à des fins de simplicités une expression continue même si celle-ci n’est pas rigoureuse dans le contexte discrétisé des valeurs obtenues par éléments finis.

Compte tenu du théorème de Cauchy égalisant les contraintes de cisaillement de facettes adjacentes, ce torseur se réduit à 8 composantes :

Considérer que cette plaque est une dalle isotrope conduit à négliger les contributions des efforts de membrane.

Il ressort donc un torseur à 5 composantes (en considérant que la notation M_x désigne le moment sollicitant les aciers dans la direction e_x) qui sont les actions mécaniques internes d’un élément plaque :



(1) Effort tranchant hors plan (direction z) pour la facette de normale x,

(2) Effort tranchant hors plan (direction z) pour la facette de normale y,

(3) Moment de flexion faisant travailler les aciers dans la direction x,

(4) Moment de flexion faisant travailler les aciers dans la direction y,

(5) Moment de torsion de la section de la dalle de normal x ou y.


3. Méthode des coupures : éléments de réduction élémentaires (EF) à éléments de réduction structurels de type poutre

Dans certains les contraintes du projet peuvent amener à modéliser intégralement une structure par plaque y compris dans des zones où celles-ci ne sont pas appropriées tels que les trumeaux et les linteaux.

Dans ces zones, il incombe de reconstituer un torseur de type poutre pour apprécier plus proprement le comportement local, en particulier pour estimer le ferraillage requis.


Une pratique répandue répondant à ce besoin est la réalisation de coupures et l’estimation des efforts sur cette coupure.

La coupure doit est choisie judicieusement de façon à adhérer aux hypothèses d’Euler-Bernoulli concernant les sections de type qui doivent rester planes et droites. C’est généralement le cas des trumeaux et linteaux.

Exemple de coupures dans un linteau ou un trumeau passant ou ne pas passant pas par des nœuds


On construit dans le repère local de la coupure les éléments de réduction de poutre à part des éléments de réduction de coque. Ceux-ci peuvent être des efforts aux nœuds ou des efforts évalués en un point quelconque de l’élément coïncidant avec la ligne de coupure et évalués au moyen des fonctions de forme élémentaires.

On écrit dans le cas de l’exemple de la figure 5 pour la ligne de coupure du linteau où le repère global et le repère de la ligne coïncident :

(1) Effort normal à la section de la poutre assimilée i.e. la coupure,

(2) Effort tranchant dans le plan,

(3) Effort tranchant hors plan,

(4) Moment de flexion dans le plan,

(5) Moment de flexion hors plan,

(6) Moment de torsion.

$translationBooks

E.2 Grandeurs en dynamique

E.2 Grandeurs en dynamique


1) Analyses temporelles

Le post-traitement des grandeurs temporelles ne présente pas de difficulté majeure tant que l’on ne s’intéresse pas à cibler un instant précis où à caractériser une valeur unique représentative de l’ensemble des instants couverts par l’analyse.

Pour les cas où le chargement est maîtrisé et ne présente pas de caractère aléatoire, une seule analyse peut suffire. S’il est privilégié d’extraire une grandeur scalaire représentative de l’ensemble des instants, le choix de la norme incombe à l’ingénieur qui analyse les résultats et qui devra être en mesure de le justifier. Une caractérisation statistique peut également être pertinente (percentile typiquement).


Pour le cas où l’ingénieur s’intéresse à extraire un ensemble de grandeurs constituants un vecteur, la problématique de la concomitance des grandeurs normées (valeurs absolue, maximale, minimale, …) se pose. L’instant où un moment est maximal en un point pour un torseur ne coïncide pas nécessairement avec l’extremum d’un effort de cisaillement.

Lorsque le chargement revêt un caractère aléatoire, il importe de multiplier le nombre de cas de calculs en intégrant le caractère aléatoire du chargement. A titre d’exemple, si on considère N cas de calculs et que l’on construit pour une quantité d’intérêt gi(t) la valeur de dimensionnement Gi correspondante pour un cas i (il peut s’agit de la valeur positive max et aussi de la valeur négative max pour les intégrer ensuite dans une démarche de concomitance d’actions).

On considère :

  • La valeur de dimensionnement retenue qui peut être (à titre d’exemple) :

  • Sa moyenne, 
  • Son écart type : 
  • L’estimation de la moyenne de la population de la valeur G pour les N résultats de calculs : 


Où λ(N) est calculé à partir de la variable de Student pour N échantillons (cas de calculs) pour un niveau de confiance 95% (unilatéral) tel que :


Le tableau ci-dessous fournit des estimations de la moyenne pour des nombres de cas de calculs variables.

N calculs t0.05;N-1 λ(N) N calculs t0.05;N-1 λ(N) N calculs t0.05;N-1 λ(N) N calculs t0.05;N-1 λ(N)
2 6.3137 4.46 15 1.7613 0.45 28 1.7033 0.32 41 1.6839 0.26
3 2.9200 1.69 16 1.7531 0.44 29 1.7011 0.32 42 1.6829 0.26
4 2.3534 1.18 17 1.7459 0.42 30 1.6991 0.31 43 1.6820 0.26
5 2.1318 0.95 18 1.7396 0.41 31 1.6973 0.30 44 1.6811 0.25
6 2.0150 0.82 19 1.7341 0.40 32 1.6955 0.30 45 1.6802 0.25
7 1.9432 0.73 20 1.7291 0.39 33 1.6939 0.29 46 1.6794 0.25
8 1.8946 0.67 21 1.7247 0.38 34 1.6924 0.29 47 1.6787 0.24
9 1.8595 0.62 22 1.7207 0.37 35 1.6909 0.29 48 1.6779 0.24
10 1.8331 0.58 23 1.7171 0.36 36 1.6896 0.28 49 1.6772 0.24
11 1.8125 0.55 24 1.7139 0.35 37 1.6883 0.28 50 1.6766 0.24
12 1.7959 0.52 25 1.7109 0.34 38 1.6871 0.27 51 1.6759 0.23
13 1.7823 0.49 26 1.7081 0.33 39 1.6860 0.27 55 1.6620 0.22
14 1.7709 0.47 27 1.7056 0.33 40 1.6849 0.27 60 1.6558 0.21

2) Analyses spectrales sismiques

Les analyses spectrales sismiques fournissent des grandeurs représentatives de la moyenne d’un extremum où cours du temps.

L’objet de ce paragraphe est d’attirer l’attention du lecteur sur des précautions importantes à prendre pour éviter de commettre des erreurs importantes dans le calcul de résultats issus d’un cumul quadratique.


Comme cela est évoqué dans le chapitre dynamique, il est possible de cumuler quadratiquement simplement ou de manière complète les réponses spectrales de chaque mode. Le résultat de ces cumuls est une valeur positive.

Il importe de considérer que toute opération que l’on souhaite réaliser sur une quantité d’intérêt qui fait l’objet d’un cumul quadratique (simple ou complet) des réponses spectrales de chaque mode doit être réalisée avant le cumul.


Prenons une illustration simple avec la différence entre le déplacement de deux nœuds A et B cumulés simplement sur 2 modes 1 et 2 (une illustration avec une combinaison CQC serait juste plus lourde en écriture en faisant intervenir les termes croisés) :

En toute rigueur, la différence évoquée à titre d’exemple s’écrit :

Il est donc assez évident de voir qu’estimer a posteriori du cumul la différence est largement erroné :


La première valeur est toujours positive mais elle intègre des écarts algébriques de quantités sur un même mode. La seconde expression peut en outre conduire à sous-estimer drastiquement une réponse extremum. On peut par ailleurs démontrer en ayant recours à l’inégalité de Cauchy-Schwartz que la première estimation est la borne supérieure de la valeur absolue du calcul erroné présenté dans la seconde formule.


On peut également prendre une illustration plus complexe avec l’estimation d’une contrainte de von Misès sur la base des contraintes principales pour 2 modes :


Une estimation a posteriori du cumul des contraintes de von Misès est largement erronée :


Exemple d’application spécifique :
 estimation de l’ouverture d’un joint entre deux bâtiments

En ingénierie sismique, il est nécessaire d’estimer l’ouverture d’un joint entre deux bâtiments sous séisme pour garantir que l’on se prémunira de tout risque d’entrechoquement. La pratique admise consiste à calculer indépendamment les valeurs extrêmes de déplacement de l’enveloppe de chaque bâtiment sur la base d’un cumul quadratique puis d’évaluer l’ouverture maximale du joint en calculant la différence entre ces deux valeurs positives. On prend ainsi en compte une opposition de phase défavorable des réponses maximales.

$translationBooks

E.3 Cas spécifique du béton armé

E.3 Cas spécifique du béton armé


1) Aperçu des méthodes usuelles de détermination du ferraillage pour des éléments de plaque


Il existe principalement 3 méthodologies applicables aux plaques pour lesquels le lecteur pourra consulter les références associées :


2) Exemple de la méthode de Capra-Maury

On définit un ensemble de facettes, centrées au point de calcul du code aux éléments finis. Il peut s’agir d’un nœud, d’un point de Gauss ou d’un point quelconque où les efforts sont interpolés.

Sa normale tourne de ce point dans le plan tangent au feuillet moyen. La facette est repérée par l’angle θ que fait sa normale avec l’axe OX du repère de l’élément (voir figure 2.1-a). L’angle θ est discrétisé régulièrement de -90°à +90° (ici avec un pas de 5°). Les axes Ox et Oy sont les axes des nappes d'armatures.

Facette de référence parallèle à la section de poutre équilibrée en flexion composée

Pour chacune de ces facettes, on évalue le moment de flexion (M), l’effort de membrane (N) et l'effort tranchant (V) qui s’y appliquent en fonction des tenseurs des efforts à l’aide des équations :


Par un calcul en flexion composée, on détermine les sections d'acier en nappes supérieure et inférieure AS(θ) et AI(θ) d'acier requise dans la direction θ pour équilibrer la section dans le contexte réglementaire de béton armé retenu.

Les efforts résistants dans la direction θ des deux nappes peuvent être évalués à l’aide des expressions

où fyd représente la contrainte maximale admissible de l’acier (identique dans les deux directions).

La résistance est assurée si l’effort résistant est supérieur à l’effort appliqué, ce qui s’écrit :


Ainsi, en considérant un repère orthonormé comportant AXS en abscisse et AYS en ordonnée, on a à résoudre finalement pour les ferraillages supérieur et inférieur :

 pour tous les angles θ

 pour tous les angles θ

et

 et  minimum.

Les inégalités sur la résistance définissent pour chaque valeur de θ un demi-espace limité par une droite de pente négative qui traduit un domaine de validité tel que cela est représenté figure suivante.

Domaine de résistance pour une facette θ

En parcourant toutes les valeurs de θ, on obtient le domaine de validité indiqué sur la figure suivante, délimité par la ligne brisée ABCD …

Domaine de résistance pour l’ensemble des facettes


Pour chaque point P du domaine de validité, la section totale des armatures peut être obtenue en projetant le point P en Q sur la première bissectrice. La distance OQ représente alors la valeur  avec AS = AXS + AYS.

On constate donc que l’optimum de ferraillage correspond à l'un des 36 points (compte tenu du pas de rotation des facettes retenu si on prend une facette tous les 10 degrés par exemple) de la frontière (illustration par les 4 points de la ligne brisée ABCD …) dont la projection sur la première bissectrice est la plus proche de l’origine des axes. La recherche de ce point peut être effectuée par une méthode de type « dichotomie ».

3) Méthodes des bielles et tirants à partir d’un résultat de calcul aux éléments finis

En présence d'éléments structurels soumis à des charges ponctuelles importantes ou présentant des modifications brusques de leur section et de leur géométrie, les méthodes classiques d'analyse des sections planes ne sont pas satisfaisantes. Ces lieux sont généralement ferraillés à l'aide de règles de bonnes pratiques basées sur l'expérience ou sur des directives empiriques. La méthode bielles-tirants (MBT) est une procédure de conception rationnelle pour les éléments structurels complexes ; la procédure a une base en mécanique, mais elle est assez simple pour être facilement appliquée dans la conception.

De manière générale, la MBT implique l’idéalisation d’un élément structurel complexe en une simple structure capable de représenter le cheminement des contraintes au sein de l’élément.

Le treillis est composé de bielles modélisant les champs de compression du béton, de tirants représentant une armature en acier élastique et de nœuds représentant les zones localisées où les éléments s’interconnectent ou les zones où l’acier élastique est ancré dans le béton. Les bielles et les tirants ne portent que des forces uniaxiales. Ce mécanisme doit être stable et équilibrer correctement les charges appliquées.

La ruine de la structure est dictée par la rupture d'un ou plusieurs tirants ou également définie par des contraintes de compression excessives à l'intérieur des bielles ou des nœuds. Idéalement, seul le premier mode de défaillance devrait se produire.

Exemple de région D et de modélisation du système pour obtention des champs de contraintes puis des efforts aux appuis du treillis BT


On applique cette méthode pour des régions dites D. Pour caractériser ces régions, on considère que la répartition des déformations sur la profondeur de l’élément présente un profil non linéaire.

Par conséquent, les hypothèses sous-jacentes à la procédure de conception de coupe sont invalidées. On ne peut y décliner le principe d’équivalence


Selon le principe de Saint Venant, une analyse des contraintes élastiques indique que les contraintes dues aux forces axiales et à la flexion se rapprochent d’une distribution linéaire à une distance approximativement égale à la profondeur de l’élément, h, à l’écart de la discontinuité.

En d'autres termes, une distribution de contrainte non linéaire existe dans la profondeur d'un membre à partir de l'emplacement où la discontinuité est introduite. Ensuite, on peut affirmer que les régions D sont supposées s’étendre jusqu’à une distance h de la charge appliquée et des réactions de support. En général, une région d'un élément structurel est supposée être dominée par un profil non linéaire, ou une région D, lorsque le rapport étendue / profondeur, a / h, est inférieur à 2 ou à 2,5. L'étendue de cisaillement, a, est définie comme la distance entre la charge appliquée et le support le plus proche dans les éléments simples.


L’approche à suivre pour définir un treillis bielle-tirant peut être résumé dans la figure 10.

Principe général de conception par MBT

Les méthodes de vérifications sont codifiées dans l’Eurocode 2 partie 1-1.


La démarche présentée dans la figure 10 est complexe à mettre en œuvre et dépend fortement de l’ingénieur qui la met en œuvre.

Des approches se développent de plus en plus dans ce registre pour chercher à automatiser la démarche. En France, le code aux éléments code_aster intègre un opérateur CALC_BT qui permet de rendre semi-automatique la démarche sur la base :

  • d’une analyse des champs pics locaux des champs de contraintes principales majeures et mineures,
  • d’un découpage de la région D modélisée par un pavage de Voronoï,
  • par projection des directions moyennes des contraintes principales dans les pavés de Voronoï
  • un ensemble de procédés d’optimisation.

Exemple de référence à gauche – Solution obtenue automatiquement par l’opérateur CALC_BT de code_aster

Cette méthode nécessite un niveau d’expérience important et un contrôle par un expert.

$translationBooks

Chapitre F - Calculs géotechniques

Chapitre F - Calculs géotechniques


Introduction

Les méthodes classiques de conception et de dimensionnement des ouvrages géotechniques concernent, en majorité, l’analyse de la résistance vis-à-vis de la rupture d’un ouvrage isolé. Ces méthodes analytiques ou semi-analytiques ne prennent en compte que des géométries très simples, et fournissent peu ou pas d’informations sur les déformations du terrain qui entoure les ouvrages.

L'utilisation de plus en plus intensive de l’espace et du sous-sol urbains, occupés par des ouvrages variés proches les uns des autres, impose de maîtriser les interactions entre ouvrages. Le concepteur d’un ouvrage doit justifier que les déplacements induits par sa construction demeurent en-deçà de seuils fixés par le maître d’ouvrage. Les méthodes traditionnelles ne répondent pas à ce besoin, et cela explique le recours de plus en plus fréquent à la modélisation numérique, à l’aide de logiciels dédiés à la géotechnique et adaptés au travail en bureau d’études. De manière plus précise, la modélisation numérique est utilisée dans deux situations différentes :

  • au stade du projet, pour justifier un dimensionnement, lorsque les méthodes traditionnelles sont difficiles ou impossibles à mettre en œuvre ;
  • comme outil d'expertise, pour étudier le comportement d'un ouvrage endommagé, identifier les phénomènes responsables d’une pathologie, et justifier l’utilisation d’une méthode de confortement.


Selon le cas, on utilise des modèles plus ou moins complexes, en tenant compte de l'incertitude sur le comportement des sols naturels et de leur variabilité spatiale. Cette incertitude se traduit par la difficulté de choisir un modèle de comportement et d’en déterminer les paramètres. L'utilisateur doit souvent faire un choix entre un modèle robuste, dont il comprend bien le fonctionnement, mais qui ne rend pas compte de toute la complexité du comportement du sol, et un modèle potentiellement plus fidèle au comportement réel des sols, mais comportant de nombreux paramètres dont le rôle est parfois difficile à cerner et qui sont difficiles à mesurer.

La philosophie des règles de dimensionnement encourage à la première approche, et à utiliser des modèles simples et robustes, et à assurer la sécurité du dimensionnement en affectant les résultats de facteurs appropriés et en effectuant des études paramétriques pour estimer la sensibilité des résultats vis-à-vis des paramètres, mais cette approche fait courir le risque de trop simplifier le problème rencontré, ou de conduire à des dimensionnements excessivement conservatifs, et inutilement coûteux.


L'autre approche consiste à utiliser des modèles de comportement qui tentent de mieux représenter les différents aspects du comportement des sols. Un très grand nombre de modèles ont été proposés, mais leur utilisation pratique reste difficile. D’autre part, il faut aussi lutter contre l’illusion d’avoir le modèle « universel » qui décrit tous les phénomènes. Il est nécessaire de cerner les limites du modèle retenu : un modèle élastique parfaitement plastique ou avec un écrouissage isotrope, même avec une surface de charge et une loi d'écrouissage complexes, ne permet pas de prédire une accumulation progressive de déformation dans un massif de sol.

L'utilisation des éléments finis demande donc de la part de l'utilisateur un recul critique sur le modèle de comportement du sol. Mais ce n'est pas le seul aspect important : le résultat des calculs peut dépendre de manière essentielle de la géométrie tridimensionnelle de l’ouvrage et de celle des couches de sol. Il est parfois (pas toujours) justifié de simplifier la représentation du comportement des sols et de privilégier celle du processus de construction. En tout état de cause, il conviendra de contrôler soigneusement les résultats des simulations numériques.


Enfin, il est clair que les ouvrages géotechniques sont constitués de terre, au contact de la terre, ou enterrés totalement ou partiellement, et que la modélisation numérique doit donc représenter l’interaction mécanique entre le sol et la structure.

Ce chapitre est organisé comme suit : il récapitule les principales spécificités des calculs numériques en géotechnique, avant de présenter succinctement les principes de vérification adoptés dans l’Eurocode 7. Une fois ces éléments présentés, on propose des recommandations visant à une bonne pratique des calculs par éléments finis en géotechnique. Il aborde également les spécificités des calculs en dynamique.

F.1 Aspects géométriques

F.1 Aspects géométriques

F.2 Non-linéarités matérielles

F.2 Non-linéarités matérielles

F.3 Interactions sol-structure

F.3 Interactions sol-structure

F.4 Effets hydrauliques

F.4 Effets hydrauliques

F.5 Incertitudes et recommandations

F.5 Incertitudes et recommandations

F.6 Aspects normatifs : Principes de l'Eurocode 7

F.6 Aspects normatifs : Principes de l'Eurocode 7

F.7 Modélisation en dynamique

F.7 Modélisation en dynamique

F.8 Échelles caractéristiques

F.8 Échelles caractéristiques


Conclusions et références

Ce chapitre présente un panorama des difficultés posées spécifiquement par les calculs par éléments finis en géotechnique, en statique ou en dynamique. Des précautions particulières doivent être prises pour la mise en œuvre des calculs, concernant le maillage, les conditions aux limites, le phasage, le choix des modèles de comportement et le calage des paramètres.

Moyennant une certaine habitude, et la préoccupation constante de contrôler les résultats obtenus, on est en mesure de traiter une large gamme de problèmes et d’ouvrages. Il reste cependant des situations qu’il est difficile et/ou coûteux de traiter numériquement.


En premier lieu, des difficultés peuvent provenir du fait que la structure mathématique des problèmes étudiés pose problème : en dynamique, le choix de pas de temps cohérents avec toutes les données du calcul reste délicat ; en statique, la gestion de plusieurs non-linéarités de natures différentes (contact unilatéral, endommagement) peut conduire à des difficultés numériques difficiles à surmonter.

Il faut aussi rappeler que la physique de certains phénomènes reste mal maîtrisée ou mal décrite, et le traitement numérique reflète cette situation : c’est le cas pour l’initiation et l’évolution des glissements de terrain, les conditions d’exécution des travaux comme le compactage des matériaux de remblai, l’état initial des contraintes, de la teneur en eau et d’autres paramètres liés à l’histoire des matériaux en place, etc.

Pour les références spécifiques, se reporter à la bibliographie :


Lien vers la bibliographie

$translationBooks

F.1 Aspects géométriques

F.1 Aspects géométriques

L’un des atouts de la méthode des éléments finis réside dans la possibilité de décrire la géométrie exacte des ouvrages, y compris au cours des différentes phases de construction. Des outils de pré-traitement proches de la CAO permettent de générer facilement des géométries très complexes. Une des spécificités des ouvrages géotechniques réside dans le fait qu’il faut en général prendre en compte dans le maillage le massif de terre qui constitue ou qui entoure l’ouvrage lui-même.


1) Limites du domaine d'étude

Une première difficulté consiste à cerner les limites du domaine pris en compte dans l’étude. Pour un ouvrage géotechnique, les limites horizontales et la limite inférieure du domaine d’étude sont rarement précisément déterminées : on limite l’étendue du domaine d’étude par des plans verticaux, dont la position est fixée généralement suivant des règles empiriques.

En déformation plane par exemple, la position à laquelle on fixe la limite inférieure du maillage a une influence directe sur le tassement calculé pour une semelle filante ou au-dessus d’un tunnel. Cette influence est nette dans le cas d’un massif de sol élastique linéaire homogène. On peut la réduire en prenant en compte des modules élastiques qui augmentent avec la profondeur, mais elle reste néanmoins susceptible d’induire une erreur significative dans les déplacements calculés. Le cas idéal est celui où l’on a reconnu un substratum rigide à une profondeur donnée, ce qui suppose que les reconnaissances aient été conduites jusqu’à une profondeur suffisante. Pour un tunnel, par exemple, il serait souhaitable, pour les besoins de la modélisation, de mener des reconnaissances bien au-delà de la profondeur de l’axe : ce n’est généralement pas le cas dans les projets réels.


Dans les directions latérales, la prise en compte d’un domaine trop peu étendu peut également modifier significativement la réponse du modèle numérique. Des déplacements bloqués conduisent à surestimer la raideur du massif, des conditions de type « contact lisse » conduisent à l’inverse à surestimer les déplacements. Le choix des dimensions du maillage adaptées à un ouvrage reste un problème largement ouvert, même si certains auteurs ont proposé des règles pratiques, qui ne doivent cependant pas être prises comme des prescriptions absolues (voir Mestat et Prat, 1999).

Le choix de l’étendue du domaine pris en compte dans le maillage est donc un point important de la modélisation des ouvrages géotechniques, même pour des analyses relativement simples en statique. Dans le cas des calculs dynamiques, la question de l’étendue du domaine maillé soulève des difficultés spécifiques et fait partie intégrante de la stratégie de modélisation. On détaillera ce point dans la section 8 ci-dessous.


2) Hétérogénéités du sol

Dans certaines régions, comme celles de Londres ou de Francfort, la géologie permet de considérer le sol au voisinage de l’ouvrage comme homogène (au sens où son comportement mécanique et hydraulique peut être représenté par un modèle unique). Cependant, il est fréquent, dans d’autres contextes, en particulier dans la région parisienne, que le domaine étudié comporte plusieurs couches de terrains présentant des natures et des caractéristiques (en particulier mécaniques) bien différentes. L’élaboration d’un modèle commence donc, comme pour les méthodes traditionnelles, par une étude détaillée des couches de sol dans la zone intéressée par l’ouvrage. Il ne s’agit pas de reproduire la géométrie exacte des couches géologiques (qui peuvent être de faible épaisseur localement), mais de définir des ensembles géotechniquement homogènes.


3) Discontinuités

Une particularité importante des calculs géotechniques réside dans la présence, au sein des massifs, de surfaces de fracture préexistantes à la mise en place de l’ouvrage ou du chargement étudié. Elles produisent une discontinuité du déplacement entre les blocs situés de part et d’autre de la surface de fracture. La méthode des éléments finis est plutôt adaptée à la recherche de champs de déplacement continus, et la prise en compte de ce type de discontinuité demande la mise en œuvre de techniques spécifiques (on utilise généralement des éléments particuliers), voire d’utiliser une autre méthode de calcul comme la méthode des éléments distincts.


L’analyse a posteriori de glissements de terrain, par exemple, est un exercice délicat. Il est à l’heure actuelle extrêmement difficile de prévoir l’apparition et le développement d’une surface de rupture. On en est réduit le plus souvent à prendre en compte une surface existante, dont on a pu reconnaître la position de manière plus ou moins précise à l’aide de dispositifs ad hoc (par exemple des inclinomètres qui suivent la déformation d’un versant).

Les massifs rocheux sont souvent traversés par un grand nombre de fractures présentant des orientations pratiquement parallèles à une ou deux directions privilégiées (diaclases régionales). La distribution des fractures est aléatoire, ou en tout cas impossible à caractériser complètement à l’échelle du massif, et diffuse. Si l’on s’intéresse au comportement global du massif, on peut proposer de le modéliser comme un milieu continu, en lui affectant un modèle de comportement intégrant, à l’échelle du calcul, l’effet des discontinuités : on recourt en général à des méthodes d’homogénéisation.

On peut aussi avoir à prendre en compte une discontinuité de grande ampleur dans un massif rocheux (faille), que l’on peut traiter comme les surfaces de rupture des glissements de terrain.


4) Un système matériel "ouvert" et des techniques de construction à modéliser

Comme on l’a dit dans le chapitre introductif, une spécificité des calculs de génie civil tient à la nécessité de prendre en compte des phases de construction, telle que le déblaiement ou le remblaiement, la prise d’un massif de béton, la mise en tension de câbles, etc. La prise en compte de ces phases de construction dans le cadre de la méthode des éléments finis n’est pas forcément simple et immédiate : la méthode consiste à se ramener au calcul d’une matrice de rigidité, d’un vecteur de forces nodales, et à résoudre le système obtenu en tenant compte des conditions aux limites. Pour modéliser les phases de construction, on est amené à enchaîner des calculs, en prenant en compte le changement de raideur de certains éléments, la disparition ou le changement de nature de certains appuis, des changements de point d’application des chargements, etc.


La difficulté consiste donc à proposer des techniques de simulation permettant de prendre en compte un grand nombre de dispositions constructives dans le cadre relativement étroit de la méthode des éléments finis. Il revient à l’utilisateur de discerner si les outils de modélisation proposés par les logiciels reflètent correctement les phénomènes mis en jeu.

Dans le cas des tunnels creusés au tunnelier, les sollicitations appliquées au terrain au cours des différentes phases d’avancement sont complexes, le terrain se refermant sur l’anneau constitué par les voussoirs lorsque le tunnelier avance. La mise en place de pieux par battage est un autre exemple de problème qu’il est difficile de ramener au cadre de la méthode des éléments finis.

$translationBooks

F.2 Non-linéarités matérielles

F.2 Non-linéarités matérielles

En géotechnique, il est très rare de pouvoir se limiter à un comportement linéaire pour l’étude d’un ouvrage (à l’exception de certaines analyses dynamiques). Pour autant, il peut être utile de réaliser dans un premier temps un calcul linéaire, pour vérifier que la géométrie et les conditions aux limites sont correctes et se faire une première idée de la déformation que le chargement peut provoquer. Cette première idée peut cependant être complètement fausse : dans le cas d’une excavation devant une paroi moulée par exemple, la cinématique calculée avec un comportement linéaire est nettement différente de celle que l’on observe.


1) Lois de comportement

Même en se limitant aux cas parfaitement saturé ou parfaitement sec, le comportement des sols est complexe. Dans la pratique, on utilise le plus souvent des modèles élastoplastiques, qui donnent une relation entre contraintes et déformations qui est non linéaire, mais indépendante du temps. Les effets de fluage et de viscosité peuvent être pris en compte pour des applications particulières qui le nécessitent – et si on peut accéder expérimentalement aux paramètres correspondants – mais l’utilisation des modèles de ce type reste limitée.

Parmi les modèles élastoplastiques, on recourt encore le plus souvent aux modèles élastiques linéaires parfaitement plastiques (voir l’enquête citée par Gilleron, 2016). L’utilisation de modèles élastiques non linéaires associés à un ou plusieurs mécanismes plastiques écrouissables se généralise progressivement, en particulier sous l’impulsion des éditeurs de logiciels : ils donnent des résultats nettement plus représentatifs de la réalité dans certains cas (par exemple pour l’excavation devant une paroi moulée) mais l’on ne maîtrise pas forcément très bien l’influence de chacun des paramètres de ces modèles. De manière générale, le choix d’un modèle de comportement pour les sols doit tenir compte des objectifs fixés pour le calcul, du type d’ouvrage (et du type de sollicitation auquel le sol sera soumis), du niveau de précision des reconnaissances et des essais en laboratoire disponibles.


2) Etat initial

Pour les modèles non linéaires, la raideur du matériau dépend de l’état initial des contraintes dans le massif de sol étudié (et éventuellement d’autres paramètres, d’écrouissage ou d’endommagement par exemple). La détermination des contraintes initiales a donc une influence déterminante sur les résultats. Malheureusement, les contraintes initiales sont généralement évaluées de manière très simple : elles sont soit assimilées à un champ de contrainte « géostatique » (pour un massif dont la surface est horizontale), soit obtenues en appliquant la gravité à l’ensemble du maillage à partir d’un état de contraintes nul. Ces procédés sont relativement pauvres comparés à la complexité des modèles rhéologiques utilisés pour les sols. Ils restent cependant incontournables en pratique, faute d’un meilleur moyen d’estimer les contraintes initiales dans le sol.

$translationBooks

F.3 Interactions sol-structure

F.3 Interactions sol-structure

Les ouvrages géotechniques associent souvent des couches de sol et des structures métalliques ou en béton, qui sont généralement beaucoup plus raides que le sol. L’interaction peut être limitée à quelques points d’appui de la structure sur le sol, ou être continue sur une surface significative, à l’extrados d’un tunnel, ou sur une paroi de soutènement.

L’interaction sera traitée de manière plus ou moins précise selon le cas.

 

Dans le cas des tunnels par exemple, on considère très couramment une adhérence parfaite entre le sol et la voûte d’un tunnel, en raison du mode de construction : avec les méthodes traditionnelles (séquentielles), le soutènement est constitué de béton projeté directement sur la surface de terrain découverte par l’excavation, ce qui assure en principe une bonne continuité du déplacement ; lors du creusement au tunnelier, on s’efforce d’assurer une bonne transmission des efforts entre le terrain et les voussoirs en procédant à des injections de bourrage pour remplir l’espace entre l’anneau et le terrain. Dans les tunnels anciens, les pathologies observées (ou des essais effectués dans le tunnel) peuvent cependant laisser penser que le contact est localement perdu entre la voûte et le terrain, par exemple à cause de circulations d’eau qui ont pu lessiver le terrain : la modélisation doit alors décrire de manière plus précise les conditions de contact entre le terrain et la voûte.

Les tunnels construits en tranchée couverte constituent un problème différent, dans la mesure où l’on remblaie le terrain autour de la structure. La modélisation de cette opération peut demander de prendre en compte explicitement, à l’aide d’éléments spécifiques, l’interface entre la voute et le sol.

 

Il est courant également d’introduire une modélisation explicite de l’interface entre le sol et la structure pour les soutènements, lorsqu’on remblaie derrière un mur (le phénomène de glissement du sol à l’interface avec la structure étant analogue à celui mis en jeu pour les tranchées couvertes), ou lorsqu’on excave devant une paroi moulée par exemple, le massif de sol soutenu pouvant glisser et présenter un déplacement vertical plus important que la tête de la paroi.

La question de la modélisation d’une interface entre les sols et les structures doit être examinée au cas par cas. On peut introduire des éléments de contact ou d’interface spécifiquement destinés à représenter l’interaction mécanique entre les deux, mais ces éléments introduisent de nouveaux paramètres, qui ne sont pas forcément simples à identifier (les raideurs normale et tangentielle de l’interface). Cette approche de modélisation présente un risque : les éléments d’interface tendent à contrôler le comportement de l’ouvrage et à estomper le rôle du comportement du sol, donnant l’impression que la réponse de l’ouvrage ne dépend pratiquement plus du sol.

 
Ouvrages renforcés

Dans de nombreux cas, le sol est renforcé par des inclusions présentant des caractéristiques de raideur et de résistance très élevées. Ces inclusions sont réparties de manière discrète dans le sol et très élancées : pieux, micropieux, tirants d’ancrage, armatures de murs en sol renforcé. Cette particularité géométrique pose différentes difficultés. En premier lieu, en toute rigueur, une file de pieux n’est pas équivalente à un mur continu, et l’utilisation de calculs en déformation plane n’est pas justifiée. En pratique, on est conduit à adopter, pour le mur du calcul plan, des caractéristiques mécaniques « équivalentes » à celles de la file de pieux, moyennant des hypothèses plus ou moins difficiles à justifier. Il en va de même pour les paramètres de l’interface mécanique entre le sol et les pieux / le mur. La difficulté est la même si l’on représente le mur par des éléments surfaciques ou par des éléments linéiques de type poutre.

Pour lever cette difficulté, on peut recourir à une modélisation tridimensionnelle. Mais, à cause des dimensions de la section des inclusions, il n’est pas possible de représenter dans le maillage la géométrie réelle des inclusions dès que leur nombre dépasse quelques unités : pour un mur en terre armée, avec des armatures de section 5 mm x 45 mm, à raison de 4 à 6 armatures par écaille de 0,75 m x 0,75 m, et pour un volume dont les dimensions sont de l’ordre de la dizaine de mètres, le nombre de nœuds d’un maillage qui respecterait la géométrie réelle des inclusions et qui donnerait une discrétisation acceptable dépasse les capacités de calcul actuelles. On peut donc proposer de représenter les inclusions par des éléments unidimensionnels (avec ou sans prise en compte des effets de flexion). Cette approche est critiquable au plan théorique, parce que l’introduction d’une densité linéique de force exercée par l’inclusion dans un milieu tridimensionnel n’est pas compatible avec la représentation classique des efforts intérieurs par un tenseur de contraintes. Elle est cependant utilisée, mais il faut être prudent dans l’interprétation des résultats, au moins pour ce qui concerne les contraintes au voisinage des inclusions.

 

Une solution alternative consiste à adopter des approches de type homogénéisation pour prendre en compte l’influence des inclusions sur le comportement global de l’ouvrage. Des modèles plus ou moins complexes ont été développés et implantés dans certains logiciels.

Quels que soient les choix effectués (calcul en déformation plane ou en condition tridimensionnelle, discrétisation des inclusions – par des éléments linéiques ou non – ou approche homogénéisée), il faut représenter l’interaction mécanique entre le pieu et le sol qui se produit au niveau du contact entre le sol et la paroi latérale du pieu, et aussi entre le sol et le pied du pieu. La modélisation de l’interaction mécanique au pied du pieu est particulièrement difficile à maîtriser. La modélisation d’un pieu unique par des éléments de volume, avec éventuellement des éléments d’interface avec le sol qui l’environne, donne des résultats qui dépendent du maillage et du modèle de comportement utilisé pour le sol. Il est, au minimum, nécessaire d’utiliser un modèle qui reproduit la rupture du sol en compression si l’on s’intéresse à la rupture du pieu. Les modèles de type Mohr Coulomb ou Drucker Prager, par exemple, ne sont pas adaptés dans ce contexte.

 

Dans certaines modélisations 2D ou 3D où les inclusions sont représentées par des éléments de barre ou de poutre, on associe à l’inclusion une extrémité fictive (par exemple un élément de poutre horizontal perpendiculaire au pieu), pour tenter de mieux représenter l’interaction en pied : il semble pour le moins indiqué de réaliser des études de sensibilité sur la dimension des éléments ajoutés pour vérifier la pertinence de cette approche.

Enfin, d’autres techniques de modélisation sont disponibles, qui proposent d’intégrer explicitement un modèle d’interaction pour le frottement latéral et un autre pour l’interaction de pointe via des éléments ad hoc.

Sans vouloir entrer davantage dans les détails, on attire donc l’attention de l’utilisateur sur le fait qu’il est amené à choisir entre des techniques de simulation numériques, et des modèles, qui ont une influence directe sur les résultats qu’il obtiendra.

$translationBooks

F.4 Effets hydrauliques

F.4 Effets hydrauliques


1) Couplage hydromécanique

Une autre particularité des calculs géotechniques concerne le rôle de l’eau présente dans les sols. Lorsqu’on applique rapidement un chargement mécanique sur une couche de sol saturée, il se produit une déformation instantanée et une mise en pression du fluide au voisinage de la charge appliquée. En fonction des conditions aux limites hydrauliques, le gradient de la charge hydraulique provoque une mise en mouvement du fluide, qui conduit à une redistribution de la pression et à déformation différée du sol.

On a donc à traiter un problème de couplage hydromécanique. Un cadre théorique solide a été établi par Biot (1941) et développé par Coussy (1991). Sur le plan de la résolution numérique, le problème couplé est nettement plus délicat à traiter qu’un problème classique, pour plusieurs raisons :

  • le problème comporte, en plus des déplacements, un nouveau champ inconnu, le champ de pression de l’eau,
  • il faut préciser des conditions aux limites spécifiques au problème hydraulique (définir les parties du contour qui sont imperméables et celles où la pression est imposée),
  • la solution (en déplacement et en pression) dépend du temps,
  • la nature mathématique du problème à résoudre est différente,
  • il faut décrire de manière quantitative la perméabilité des différentes couches de sol.

Le traitement complet du couplage hydromécanique est rarement effectué. On tente généralement de se limiter à une approche découplée du problème, dans laquelle on calcule l’évolution du champ de pression en négligeant les déformations du solide, mais ce découplage a pour conséquence de sous-estimer fortement la durée sur laquelle la redistribution de pression se produit.

 
2) Sols non saturés

Le traitement des sols non saturés complique encore le problème, la transition entre les zones quasi-saturées et les zones non saturées introduisant des inconnues, des non linéarités et des paramètres supplémentaires. A nouveau, le traitement complet des sols non saturés reste rare. On préfère proposer des solutions simplifiées, en négligeant les zones partiellement saturées par exemple.

La question de l’état initial dans le cas des sols non saturés est extrêmement délicate, faute de moyens expérimentaux permettant de le caractériser in situ.

$translationBooks

F.5 Incertitudes et recommandations

F.5 Incertitudes et recommandations


1) Incertitudes

Les sources d’incertitude sont nombreuses en géotechnique. La première réside dans la relative méconnaissance de la géométrie des couches de sol constituant le massif étudié. Outre la connaissance de la géologie de la région où est implanté l’ouvrage étudié, les principales sources d’information sont les sondages effectués sur le site du projet. L’extrapolation entre les sondages, surtout s’ils sont éloignés les uns des autres et pas localisés précisément à l’emplacement de l’implantation définitive de l’ouvrage (qui peut être modifiée après les reconnaissances, par exemple dans le cas des tunnels), ne donne pas nécessairement une représentation exacte des variations locales d’épaisseurs des couches.

En zone urbaine, la présence d’hétérogénéités (caves, puits, fondations d’ouvrages antérieurs) est souvent difficile à détecter.

L’autre source d’incertitude, déjà évoquée, concerne l’état initial des contraintes (et des pressions interstitielles éventuellement) dans le massif de sol. Elle peut avoir une influence majeure sur les résultats du calcul : c’est particulièrement clair dans le cas des tunnels, où les chargements pris en compte dépendent des contraintes initiales.

Enfin, le choix des modèles de comportement et la détermination des paramètres de ces modèles introduisent une incertitude importante sur la représentativité des calculs : si le modèle de comportement ne capte pas un phénomène qui contrôle le comportement de l’ouvrage, le résultat peut être qualitativement et quantitativement très éloigné de la réalité.

 

2) Recommandations

De manière générale, l’utilisateur doit être conscient des objectifs du calcul qu’il entreprend : la démarche est différente selon qu’on cherche à justifier un dimensionnement ou à évaluer l’influence de certaines dispositions constructives (le nombre et la position des butons par exemple).

Il faut également être conscient des choix de modélisation sur lesquels le calcul repose (même si ces choix sont parfois en partie imposés par le logiciel qu’on utilise). Il faut en particulier être en mesure d’identifier les phénomènes à prendre en compte, ce qui conduit à choisir une analyse quasi-statique ou dynamique, avec ou sans prise en compte du couplage hydromécanique, etc.

Il convient de choisir entre un calcul 2D et un calcul 3D. Les calculs tridimensionnels restent pour le moment rares, à cause du temps de préparation des calculs. Pour certains problèmes, cependant, il ne fait pas de doute que des calculs bidimensionnels ne peuvent donner qu’une indication très pauvre du comportement de l’ouvrage étudié, quel que soit le soin apporté à la détermination des paramètres de sol. Ainsi par exemple, l’étude de la stabilité du front de taille d’un tunnel ne peut pas vraiment être envisagée en dehors d’un contexte tridimensionnel. Il en va de même pour l’étude du renforcement du front de taille d’un tunnel par des boulons. Le développement de pré-processeurs dédiés à des applications particulières devrait permettre de recourir plus facilement à des calculs tridimensionnels et améliorer la représentativité de nombreuses analyses par éléments finis.

 

En géotechnique, le choix des paramètres de sol doit faire l’objet d’une attention particulière : il pourrait faire l’objet d’un livre entier. La plupart des modèles de comportement évolués ne s’accompagnent pas d’une procédure d’identification des paramètres détaillée et robuste, principalement parce qu’on ne sait pas résoudre les équations du modèle même pour un problème simple comme celui de la compression à l’appareil triaxial. On est donc amené à caler les paramètres en faisant en sorte que la modélisation d’essais triaxiaux par exemple donne des résultats en accord satisfaisant avec les résultats d’essais qu’on peut avoir. On procède par tâtonnements, et l’accord obtenu étant évalué de manière subjective (parce qu’on peut choisir de mieux reproduire une partie ou une autre des courbes expérimentales), il n’est pas garanti que deux utilisateurs retiennent les mêmes valeurs des paramètres à partir des mêmes essais. De ce point de vue, tous les modèles de comportement n’ont pas les mêmes qualités : certains disposent d’un grand nombre de paramètres qui ont chacun une influence sur un aspect particulier de la réponse du sol (mais qui n’apparait pas forcément dans les résultats d’essais dont on dispose) ; d’autres modèles ont au contraire un nombre relativement réduit de paramètres, mais chacun d’eux peut modifier simultanément plusieurs aspects de la déformabilité du sol : cela ne facilite pas le calage.

La dernière recommandation que l’on doit absolument garder à l’esprit est que l’utilisateur doit vérifier autant que possible son calcul. Il n’existe pas encore d’outils généraux pour mesurer la qualité d’un calcul : des travaux de recherche visent à fournir des estimateurs d’erreur a posteriori mais leur utilisation en géotechnique reste rare. Il est donc nécessaire de regarder en détail les résultats du calcul : certaines incohérences sont parfois faciles à détecter. En cas de doute, il est utile de faire contrôler ses résultats par un œil extérieur. En tout état de cause, il est vivement indiqué de réaliser des études paramétriques pour se faire une idée de l’influence de certains facteurs, en particulier des paramètres de sol, si l’on n’est pas certain que leur influence sur les résultats est modérée et qu’on a pu déterminer leur valeur avec une précision acceptable.

$translationBooks

F.6 Aspects normatifs : Principes de l'Eurocode 7

F.6 Aspects normatifs : Principes de l'Eurocode 7

La modélisation numérique a toujours eu un lien particulier avec les normes de calcul qui se focalisent essentiellement sur la vérification d’états limites ultimes avec le calcul d’un coefficient de sécurité ou d’un équilibre de forces intégrant des coefficients partiels. En effet, la modélisation numérique fournit avant toute chose des valeurs de déplacements et de déformations et est donc un outil très adapté à la vérification des états limites de service (ELS) pour lesquels les coefficients partiels sont égaux à 1,0. Son utilisation pour la vérification des états limites ultimes (ELU) a donc paru limitée dans un premier temps.

Désormais, notamment avec les procédures de réduction des propriétés de cisaillement des terrains (par exemple, la procédure de type « c-phi reduction »), il est aisé de calculer un coefficient de sécurité. Il est aussi possible, à travers les procédures suggérées par certaines normes de calcul, notamment l’Eurocode 7, de considérer les résultats d’une modélisation numérique tant pour la vérification des états limites ultimes que pour celle des états limites de service.

L’Eurocode 7 dans sa version actuelle n’est toutefois pas forcément très clair quant à l’analyse et l’exploitation des résultats issus d’une modélisation numérique. En effet, les trois approches de calcul proposées par l’Eurocode 7 qui permettent d’appliquer des coefficients partiels sur les actions, les effets des actions (moments fléchissants et efforts tranchants dans un écran de soutènement, effort axial dans un pieu, etc.), les propriétés des terrains (c et  ou cU) et les résistances géotechniques ont été pensées pour être utilisées avec des méthodes d’équilibre limite.

Les vérifications liées à l’application des Eurocodes sont principalement transcrites sous la forme de comparaisons entre des actions ou des effets des actions et des résistances. Néanmoins, un certain nombre de publications (Potts et Zdravkovic, 2012, Tschuchnigg et al., 2015, etc.) et les compte rendus de différents groupes de travail réunis en vue de la mise au point de la seconde génération des Eurocodes permettent d’esquisser des procédures de vérifications. En premier lieu, il est important de souligner que la pondération « à la source » des propriétés des sols ou des roches, c'est-à-dire la possibilité de réduire la cohésion et l’angle de frottement avant de réaliser le calcul n’est pas admise car elle conduit à générer des résultats qui ne peuvent pas être interprétés. Différentes procédures sont toutefois utilisables et sont synthétisés dans le tableau ci-dessous.

Synthèse des différentes types de vérifications aux ELU

Type de vérifications aux ELU 1 – ELU structurels 2 – ELU géotechniques (et structurels) 3 – ELU géotechniques 4 - ELU structurels (et géotechniques)
Type de procédures Multiplication des effets des actions par 1,35 Réduction des propriétés de cisaillement des terrains Estimation de la résistance mobilisée autour d’un ouvrage géotechnique spécifique (pieu, tirant, etc.) Augmentation des charges appliquées sur l’ouvrage géotechnique
Commentaires A combiner avec la vérification des ELU géotechniques A combiner avec la vérification des ELU structurels A combiner avec la vérification des ELU structurels A priori, cette approche se suffit à elle-même


L’approche qui tend à s’imposer consiste à réaliser un calcul avec des valeurs caractéristiques tant pour les charges que pour les propriétés des terrains pour atteindre un premier état d’équilibre. Par rapport à cet état d’équilibre, deux types de vérifications sont à réaliser. La première de type 1 est relative aux ELU structurels et consiste à multiplier les effets des actions calculés par 1,35 (c'est-à-dire le coefficient habituellement considéré dans les Eurocodes pour les charges permanentes défavorables). La seconde de type 2 est relative aux ELU géotechniques et consiste à réduire de manière progressive les propriétés de résistance au cisaillement des terrains modélisés pour mettre en évidence un mécanisme de rupture. Le facteur de réduction appliqué peut être considéré comme un coefficient de sécurité mais son interprétation peut être sujette à discussion. En effet, lors de la réduction de propriétés de cisaillement, différents éléments peuvent interagir comme les fonctions d’écrouissage ou les règles d’écoulement. Dans le cas de calculs associant des éléments volumiques et des éléments de structure de type « barre » ou « poutre », cette vérification peut aussi être utilisée pour vérifier les ELU structurels. En effet, les déplacements accumulés au cours de la procédure de réduction des propriétés de résistance au cisaillement engendrent des efforts dans les éléments de structure dont l’interprétation n’est pas complètement partagée par l’ensemble des chercheurs ou ingénieurs. L’utilisateur doit donc être très prudent sur la manière donc ce type de procédures est implémenté. La figure 7.1 présente comment s’enchaînent les vérifications de type 1 et 3 dans le cas d’un calcul phasé.


Les ELU géotechniques peuvent être aussi appréhendés de manière spécifique en considérant la résistance mobilisée autour de l’ouvrage géotechnique à dimensionner, un pieu, un tirant d’ancrage, etc. et à la comparer à la résistance mobilisable fournie par l’application des normes, par exemple, les règles pressiométriques pour le calcul de la portance des pieux : c’est la procédure de type 3 et elle doit forcément être associée à la procédure de type 1 en ce qui concerne les ELU structurels.

Dans certaines configurations, notamment des pieux ou des semelles, il est aussi possible de faire croître les actions appliquées sur l’ouvrage géotechnique à dimensionner pour directement obtenir l’effort maximal applicable et en déduire un coefficient de sécurité : c’est l’objet de l’approche de type 4.

De manière plus générale, les méthodes permettant de justifier les ELU géotechniques doivent permettre, soit de déterminer un mécanisme de rupture et de vérifier que l’on dispose de suffisamment de marge vis-à-vis du déclenchement de ce mécanisme (c’est l’objet des vérifications de type 2 et 4), soit de comparer la résistance mobilisée à une résistance mobilisable que l’on peut calculer par ailleurs (c’est l’objet de la vérification de type 3).

Enchaînement des vérifications aux ELU dans un calcul phasé

Quand le mécanisme de rupture est associé à l’action de l’eau alors des précautions particulières sont à prendre. En effet, les procédures décrites précédemment ne permettent pas correctement de calculer un coefficient de sécurité si le mécanisme de rupture résulte exclusivement de l’action de l’eau. En premier lieu, il faut définir si le mécanisme de rupture est lié à un défaut d’équilibre mécanique (par exemple, un niveau trop haut en arrière d’un écran de soutènement) ou à un problème hydraulique (par exemple, un gradient hydraulique trop élevé).

Dans le cas d’un défaut d’équilibre mécanique, il est nécessaire de réaliser une étude paramétrique relative au niveau de la nappe de manière à pouvoir estimer son effet en mettant en œuvre les vérifications aux ELU présentées dans le tableau précédent. Dans le cas d’un problème hydraulique, il est possible de comparer le gradient calculé à un gradient critique calculé par ailleurs et aussi de vérifier que les contraintes effectives restent positives en tout point du massif (la condition sur le gradient critique étant plus conservatrice que la condition sur la contrainte effective).

$translationBooks

F.7 Modélisation en dynamique

F.7 Modélisation en dynamique

A cause de la nature continue et non-confinée des géomatériaux, le traitement des problèmes dynamiques dans le milieu « sol » diffère essentiellement du traitement des problèmes dynamiques associés aux structures classiques du génie civil. Alors que pour la plupart de ces dernières une modélisation mécanique, par le biais d’un assemblage de masses et de sources de rigidité discrètes, suffit a priori pour appréhender leur comportement dynamique, les géomatériaux doivent être traités en tant que milieux continus et non-bornés et leur réponse dynamique doit être étudiée dans le contexte des problèmes mécaniques de la propagation d’ondes.

Concernant les méthodes numériques employées en dynamiques des sols, il convient d’abord de noter que la modélisation de la propagation d’ondes nécessite de faire appel à une palette de méthodes de modélisation et de calcul, qui diffèrent des calculs géotechniques conventionnels par la méthode des éléments finis. Bien que cette dernière occupe une place importante dans l’inventaire de méthodes à la disposition de l’ingénieur pour les problèmes de dynamique des sols, elle est souvent couplée ou remplacée par d’autres méthodes de calcul numérique, mieux adaptées pour modéliser la propagation d’ondes dans des milieux non-bornés.

 

En outre, les équations qui décrivent les problèmes dynamiques font intervenir une nouvelle variable (en dehors des variables spatiales) : celle du temps. Les schémas de résolution sont formulés dans le domaine temporel, et alternativement, dans le domaine fréquentiel. Les méthodes pour les problèmes en dynamique des sols sont implémentées dans des logiciels dédiés, nécessitant à la fois une très bonne maîtrise des notions théoriques de base et aussi un savoir-faire ciblé pour ce qui est de la définition des paramètres, du calage des modèles numériques et des lois constitutives, de l’optimisation des calculs etc.

Parmi les problèmes courants de la dynamique des sols, on peut citer :

$translationBooks

F.8 Échelles caractéristiques

F.8 Échelles caractéristiques

Pour cadrer la discussion concernant les problèmes de la dynamique de sols, il faut d’abord faire la distinction parmi les échelles caractéristiques relatives à la définition des phénomènes physiques étudiés dans ces problèmes. Ces échelles sont :

  • l’échelle géologique,
  • l’échelle de site,
  • l’échelle géotechnique,
  • l’échelle de la structure génie civil.

Chaque problème en dynamique de sols fait intervenir une ou plusieurs de ces échelles, ce qui constitue déjà une première difficulté, car il faut associer dans un même modèle, des domaines dont les dimensions sont très différentes.

 

Le passage d’une échelle à une autre correspond à un type de problème spécifique, nécessitant des approches de résolution qui lui sont propres. Ainsi, dans le contexte des problèmes sismiques par exemple, le passage de l’échelle géologique à l’échelle de site donne lieu au problème de la détermination des effets locaux de site. Le passage à l’échelle géotechnique donne lieu au problème de la caractérisation du mouvement sismique en un point de contrôle (qui peut être localisé en surface du sol ou à une profondeur caractéristique). Enfin, le passage à l’échelle de la structure donne lieu au problème de l’interaction sol-structure, c’est-à-dire au problème de déterminer comment la réponse dynamique du sol est modifiée par la présence de la structure et vice versa.

Selon le problème étudié, il est également possible que l’échelle de la structure interagisse directement avec l’échelle de site ou l’échelle géologique. Ceci est surtout vrai pour les structures de grandes dimensions, telles que les ponts (échelle de structure ~ échelle de site) ou les tunnels (échelle de structure ~ échelle géologique).

Dans la suite, il convient de limiter la discussion sur les échelles géotechnique et de la structure (échelles caractéristiques pour la majorité des structures du génie civil) et donc se concentrer aux problèmes d’Interaction Dynamique Sol-Structure (IDSS).

 
1) Classification des méthodes pour l’IDSS

Pour proposer une classification des méthodes pour l’IDSS, on notera d’abord que les systèmes dont le comportement dynamique est étudié, sont constitués des deux domaines distincts :

  1. le domaine de la structure, comportant l’ouvrage étudié et éventuellement une zone de sol au voisinage de cet ouvrage : ce domaine met en évidence des dimensions bornées. Il est couramment appelé domaine de structure généralisée ou superstructure ou tout simplement structure ;
  2. le domaine de sol, entourant cette structure, dont les dimensions sont non-bornées.

Les deux domaines sont séparés par une frontière qui est appelée horizon d’interaction. Cette frontière est bien évidemment fictive et sa définition est, dans un certain sens, arbitraire. Par exemple, l’horizon d’interaction peut être identifié avec la vraie interface entre la structure et le sol, ou être placée très loin de la structure.

 

Par ailleurs, puisque le problème d’IDSS résulte en un problème de propagation d’ondes, la formulation adoptée doit satisfaire une condition de radiation. Par exemple, dans le cas d’un chargement appliqué sur la structure (c’est-à-dire que la structure est identifiée comme la source des ondes), la condition de radiation exige qu’il ne puisse pas y avoir d’ondes convergeant vers la source.

Les méthodes pour l’IDSS sont divisées en trois grandes familles :

  • méthodes de sous-structures,
  • méthodes directes,
  • méthodes hybrides.

 

Les méthodes de sous-structures sont les plus courantes dans les études d’IDSS en Génie Civil. L’horizon d’interaction dans ce cas est placé au niveau de l’interface physique entre le sol et la structure. Cette séparation permet de modéliser le domaine « structure » par des éléments finis classiques et le domaine « sol » par d’autres approches, mieux adaptées pour le traitement de domaines non-bornés, en particulier par la méthode des éléments de frontière (BEM ou une de ses variantes). La condition de radiation est formulée de manière rigoureuse à l’infini. La méthode de sous-structures acquiert son intérêt pratique, si elle est combinée avec l’hypothèse de linéarité de comportement pour les domaines « structure » et « sol ». Dans ce cas, il peut être montré, via le théorème de superposition de Kausel, que le problème global d’interaction sol-structure peut être décomposé en trois sous-problèmes, comme présenté sur la figure ci-dessous.

Théorème de superposition de Kausel dans la méthode de sous-structures

L’interaction cinématique concerne la diffraction des ondes dans le sol par la présence de la fondation et le chargement conséquent de cette dernière. L’interaction inertielle concerne l’échange d’énergie entre le sol et la structure qui peut être caractérisé en termes d’impédance. Ces deux premiers problèmes fournissent respectivement, le mouvement d’entrainement et les conditions aux limites qui peuvent être utilisés, dans le troisième sous-problème, pour l’analyse dynamique du domaine borné. Par conséquent, dans la méthode des sous-structures, l’utilisation des éléments finis n’intervient essentiellement que dans la modélisation de la structure ; on peut donc se référer directement aux chapitres concernant l’utilisation des éléments finis dans les études conventionnelles en dynamique des structures.

Un des inconvénients, quant à l’utilisation du théorème de superposition de Kausel, est lié à la nature des géomatériaux, qui présentent d’emblée un comportement mécanique (fortement) non-linéaire. Ainsi, pour pouvoir faire appel à l’hypothèse de linéarité de comportement pour le domaine « sol », il est nécessaire de procéder à un calage des paramètres viscoélastiques des géomatériaux, en conformité avec les niveaux des déformations qui y sont anticipés. La définition des paramètres viscoélastiques effectifs des couches de sol constitue la première étape pour la mise en œuvre de la méthode de sous-structures et permet l’utilisation du théorème de superposition dans des systèmes où le comportement non-linéaire reste modéré.

 

Les méthodes directes constituent la deuxième grande famille d’approches pour l’IDSS et sont surtout applicables pour les systèmes qui présentent de fortes non-linéarités de comportement. L’horizon d’interaction dans ce cas est placé très loin de la structure, englobant une grosse zone du domaine « sol », comme présenté sur la figure ci-dessous.

Principe général des méthodes directes

L’équation d’équilibre dynamique est formulée pour l’ensemble « structure + sol » et fait intervenir un vecteur de chargement extérieur qf, défini le long de l’horizon d’interaction. L’ensemble « structure + sol » est modélisé par la méthode en éléments finis et les déplacements dans les domaines « sol » et « structure » sont calculés simultanément. La résolution est effectuée dans le domaine temporel. Il est alors possible d’incorporer dans le modèle toutes les caractéristiques géométriques du problème, les hétérogénéités matérielles du sol ou de la superstructure et d’introduire les lois de comportement nécessaires pour la description des non-linéarités du système. En dépit de ces importants avantages, les méthodes directes mettent en évidence plusieurs points délicats :

  • la définition de la position optimale pour l’horizon d’interaction et la détermination du chargement extérieur à y introduire (surtout pour les problèmes sismiques) sont des tâches non triviales, nécessitant souvent un travail préparatoire conséquent ;
  • à l’opposé des méthodes des sous-structures, où la condition de radiation est formulée de manière rigoureuse à l’infini, dans les méthodes directes, la condition de radiation est formulée de manière approchée au niveau de l’horizon d’interaction. Cette hypothèse permet de s’affranchir de la modélisation du domaine non-borné, mais crée le besoin d’utilisation d’éléments spéciaux le long de l’horizon d’interaction (qui est identifié avec la frontière du domaine modélisé) ; ces éléments sont censés préserver le caractère non-borné du domaine modélisé en garantissant qu’il n’y aura pas de réflexions d’ondes parasites le long de la frontière du domaine modélisé ;
  • les discrétisations spatiale et temporelle (taille de mailles, pas de temps) dépendent des paramètres mécaniques du système modélisé et souvent la taille de mailles ou le pas de temps qui sont nécessaire pour un problème, peuvent s’avérer prohibitifs vis-à-vis du coût de calcul ;
  • le schéma d’intégration temporelle doit satisfaire les exigences de précision et de stabilité de la solution, tout en restant optimisé par rapport au coût calculatoire ;
  • les lois de comportement adoptées pour décrire le comportement non-linéaire des sols nécessitent un calage pour pouvoir bien représenter le comportement cyclique des sols ;
  • enfin, les méthodes directes restent encore extrêmement coûteuses pour des problèmes 3D. Il est donc nécessaire souvent de pouvoir faire le passage du comportement 3D à une modélisation 2D tout en préservant les caractéristiques dynamiques (masse / rigidité / caractéristiques modales) du système étudié.

 

Enfin, les méthodes hybrides forment la troisième grande famille de méthodes pour le traitement de l’IDSS. Ces méthodes sont situées entre les méthodes directes et les méthodes de sous-structures ; l’idée principale consiste à séparer le sol dans deux domaines distincts : le premier est un champ proche de la structure où l’on suppose que toutes les non-linéarités pertinentes pour le problème de l’IDSS sont développées. Le deuxième est le champ lointain où le comportement du sol n’est pas affecté par l’interaction avec la structure. Ainsi, le champ lointain peut être traité par les techniques adaptées pour les problèmes linéaires dans des domaines non-bornés (typiquement BEM avec une résolution dans le domaine fréquentiel), alors que le champ proche est incorporé dans le modèle de la superstructure et peut être traité par une méthode directe (modélisation par éléments finis – modélisation des non-linéarités – résolution dans le domaine temporel). Le principe général des méthodes hybrides est présenté sur la figure ci-dessous.

Principe général des méthodes hybrides

 

Le point délicat des méthodes hybrides est la définition de la frontière entre le champ proche et le champ lointain (horizon d’interaction) ainsi que la formulation d’un schéma qui permet de restituer de manière rigoureuse mais efficace le couplage entre les deux champs. Au plan numérique, cette opération résulte typiquement en un schéma de couplage FEM-BEM, où la composante BEM alimente la composante FEM avec les historiques des forces d’interaction et les fonctions d’impédances à déployer le long de l’horizon d’interaction.

Un cas particulier de modélisation appartenant aux méthodes hybrides concerne le développement de macroéléments pour la modélisation du comportement dynamique non-linéaire de différents types de fondations. Les macroéléments sont des éléments de liaison, qui sont situés au niveau de la fondation et qui sont dotés d’une loi de comportement permettant de décrire les mécanismes non-linéaires développés le long des interfaces sol-structure, tels que le glissement, le décollement et la plastification du sol.

Les points précédents laissent entendre qu’au sein des approches pour l’IDSS, la méthode aux éléments finis est surtout utilisée dans la famille des méthodes directes et aussi dans la famille des méthodes hybrides pour la modélisation du champ proche.

Dans les paragraphes suivants, on donne quelques indications sur la méthode de sous-structuration et on discute quelques points délicats quant à l’implémentation des méthodes directes pour les problèmes d’IDSS.

 
2) Compléments sur la méthode de sous-structuration

La première étape de la méthode de sous-structuration est le calage des paramètres viscoélastiques des géomatériaux en fonction du niveau des déformations. Le comportement non-linéaire hystérétique d’un sol sous sollicitations cycliques peut être approximé par un module de cisaillement et un amortissement visqueux équivalents. D’abord il faut définir les courbes de la réduction du module de cisaillement et de l’augmentation du coefficient d’amortissement en fonction de la déformation sur la base de résultats d’essais de laboratoire dédiés au projet (partie des reconnaissances géotechniques) ou, en l’absence d’une campagne géotechnique dédiée, sur la base de références scientifiques.

Puis, une analyse unidimensionnelle linéaire de la réponse d’une colonne de sol est à conduire à l’aide des fonctions de transfert qui, à partir de l’accélération au rocher et d’un premier jeu des propriétés des couches de sol, donnent les déplacements, les vitesses, les accélérations, les contraintes et les déformations en cisaillement dans chaque couche. Enfin, une procédure itérative est nécessaire afin d’assurer que les propriétés utilisées dans l’analyse sont compatibles avec les déformations maximales calculées dans les couches.

L’analyse inverse, soit le calcul du mouvement au rocher à partir de l’accélération en surface, est possible à travers les mêmes fonctions de transfert et prend le nom de déconvolution. La procédure itérative ne garantit pas toujours l’unicité de la solution dans la recherche des propriétés équivalentes du sol en cas des larges déformations. Ces analyses peuvent être menées avec des logiciels dédiés comme SHAKE 91 par exemple.

Les problèmes de diffraction et d’impédance peuvent être résolus dans le domaine fréquentiel par des outils de calculs adaptés, tel que les codes SASSI2010 ou MISS3D. Des méthodes simplifiées valables sous certaines conditions peuvent être utilisées, en alternative.

L’étude de l’interaction cinématique n’est d’intérêt que dans le cas de fondations profondes et est négligée dans le cas de fondations superficielles. Pour une analyse sismique, la méthode simplifiée consiste en l’imposition des déformations de sol en champ libre aux pieux de fondation. Cette sollicitation, d’origine cinématique, doit être additionnée aux sollicitations provenant de l’interaction inertielle pour le dimensionnement des pieux.

En utilisant un logiciel dédié (eg. SASSI2010, MISS3D), l’étude de l’interaction inertielle conduit au calcul de la raideur et de l’amortissement au point de contrôle de la fondation en fonction de la fréquence. Pour le calcul des fonctions d’impédance, la fondation est par défaut considérée comme une sous-structure rigide. Pour une fondation superficielle, les termes d’impédance relatifs aux déplacements et aux rotations sont retenus, tandis que pour une fondation profonde il faut prendre en compte aussi les termes de couplage entre les raideurs en déplacements et en rotations. Les valeurs des impédances, retenues pour le problème structurel d’analyse dynamique, sont calées sur la base des modes propres d’interaction sol-structure avec une procédure itérative expliquée ci-après. Une première analyse modale est conduite avec un jeu des valeurs d’impédance d’essai. Le jeu des valeurs d’impédance suivant est calculé à partir des fréquences des modes propres, identifiés comme ceux d’interaction sol-structure pour chaque degré de liberté de la fondation (supposée rigide). La deuxième analyse modale permet de corriger les fréquences des modes propres retenus et ainsi de suite. La procédure s’arrête lorsqu’on atteint la convergence des fréquences propres.

Les méthodes simplifiées pour les fondations superficielles, en général, donnent des formules pour le calcul direct de l’impédance. Les méthodes les plus connues sont celles de Newmark-Rosenblueth, de Deleuze, de l’Eurocode 7, du guide SETRA, de Veletsos et de Gazetas.

La condition principale pour la validité des méthodes simplifiées est l’hypothèse de sol homogène, monocouche, à l’exception de la méthode de Gazetas permettant la prise en compte d’un sol double-couche pour les cas des fondations de forme approximable à la forme circulaire. La méthode de Deleuze tient compte de la fréquence du mode de vibration fondamental de l’ouvrage. Les méthodes de Veletsos et Gazeats considèrent l’encastrement de la fondation dans le sol.

Pour les fondations profondes, la méthode de Gazetas et celle de l’Eurocode 8-5 donnent des formules applicables dans les trois conditions de sol suivantes :

  • sol dont le module de Young varie linéairement avec la profondeur ;
  • sol dont le module de Young varie avec la racine carrée de la profondeur ;
  • sol dont le module de Young reste constant avec la profondeur.

Une autre méthode simplifiée pour les fondations profondes est celle de Winkler nécessitant la modélisation d’une poutre reposant sur un système de ressorts et d’amortisseurs, calés sur la base des propriétés de chaque couche et représentant les pressions frontales et les efforts de frottement.

Suite à la résolution du problème d’analyse dynamique, la méthode de Winkler est utilisée aussi pour le calcul des sollicitations le long des pieux pour leur dimensionnement. Cette méthode permet de modéliser un comportement élasto-plastique pour chaque couche de sol et l’intégration d’une analyse simplifiée de l’interaction cinématique par application des déplacements de sol imposés au pieu.

 
3) Définition de la taille des mailles et du pas de temps

Le choix de la taille de mailles constitue une étape très importante dans la modélisation en dynamique des sols, puisque la finesse de discrétisation conditionne les caractéristiques des ondes qui peuvent être modélisées avec un maillage donné. Le critère habituel pour pouvoir modéliser avec une précision adéquate une onde se propageant dans le milieu sol, est de posséder dans le maillage, 10 mailles par longueur d’onde. Si on désigne par h la taille de mailles et par λ, la longueur d’onde, le critère à imposer sur la finesse de discrétisation s’écrit comme suit :

Ce critère peut encore être formulé en fonction de la vitesse de propagation minimale Vmin à prendre en compte dans le modèle et de la fréquence caractéristique maximale fmax des ondes à représenter dans le maillage. Le critère pour la taille de mailles s’écrit alors :

Pour les problèmes sismiques, la quantité Vmin correspond habituellement à la valeur minimale des vitesses effectives (compatibles avec le niveau sismique) de propagation d’ondes S pour les couches de sol modélisées. La quantité fmax correspond à la fréquence de coupure du spectre de dimensionnement (~30 Hz) ou à une fréquence maximale d’intérêt pour l’interaction sol-structure (~15 Hz).

En dehors de la taille de mailles h, le choix d’un pas de temps Δt est également conditionné par les caractéristiques mécaniques du milieu modélisé. En particulier, par analogie avec la longueur d’onde, il est nécessaire de pouvoir décrire la plus faible période propre à prendre en compte dans les analyses par au moins 10 incréments de temps (le critère de précision est identique dans les discrétisations spatiale et temporelle). Le critère sur le pas de temps peut être formulé comme suit :

où Vmax est la vitesse correspondant à valeur maximale des vitesses (effectives) de propagation des ondes S et h correspond à la taille de mailles.

A titre indicatif, dans une application avec Vmin = 150 m/s, Vmax = 500 m/s et fmax = 15 Hz, on aurait obtenu h ≤ 1 m et Δt ≤ 0,002 s. Une attention particulière doit être portée aux problèmes dynamiques non-linéaires, car les caractéristiques d’intérêt pour la définition de la taille de maille (en particulier V_min) évoluent avec le temps (assouplissement du sol en passant en régime élastoplastique).

 

4) Comportement cyclique des sols - Amortissement

On a vu que le grand intérêt pour la mise en œuvre des méthodes directes pour l’IDSS concerne la possibilité de modéliser rigoureusement le comportement non-linéaire des géomatériaux, au moins à l’intérieur du domaine englobé par l’horizon d’interaction, qui est modélisé en éléments finis. Le concepteur peut modéliser avec précision la géométrie des couches de sol et affecter à chaque formation, la loi de comportement la plus adaptée pour pouvoir reproduire le comportement cyclique anticipé.


Calage des modèles non-linéaires
 - La première étape dans cette procédure de modélisation concerne le calage de paramètres numériques intervenant dans la définition des lois constitutives choisies, de manière à que ces lois reproduisent avec précision le comportement cyclique des sols ; le comportement cyclique est d’emblée caractérisé, dans la majorité d’applications en Génie Civil, via la courbe de dégradation du module de cisaillement G en fonction de la distorsion γ (courbe G~γ). Il est rappelé que la courbe G~γ varie selon le type de sol (sable/argile/silt/gravier), la contrainte verticale effective, l’indice de plasticité etc. La procédure du calage des modèles non-linéaires intègre les étapes suivantes :

  • organisation du maillage faisant apparaître des zones régies par la même courbe G~γ ;
  • définition des paramètres élastiques (paramètres à faibles distorsions) pour chaque zone ;
  • adoption d’une loi constitutive adaptée pour chaque zone (surface de rupture, loi d’écrouissage, loi d’écoulement, comportement dans le plan déviatorique et dans l’axe hydrostatique etc.) ;
  • définition de la courbe G~γ caractérisant chaque zone en fonction -entre autres- de la contrainte moyenne effective dans la zone ;
  • calage des paramètres numériques caractérisant la loi de comportement non-linéaire de manière à reproduire la courbe G~γ de chaque zone.

La procédure du calage résulte en un tableau fournissant pour chaque zone, le type de loi constitutive adoptée et l’ensemble des paramètres numériques, nécessaires pour sa définition.


Comportement drainé et non-drainé
 - Dans le cas de la présence d’une nappe phréatique dans le domaine de sol, il est important de pouvoir définir si la réponse des géomatériaux lors de la sollicitation dynamique sera régie par un comportement en conditions drainées ou non-drainées. Cette définition dépend de la nature des sols (matériaux cohésifs ; matériaux pulvérulents) mais également de la durée du phénomène étudié. Pour les phénomènes dont la durée est limitée (de l’ordre de quelques secondes, y compris les séismes) il est souvent difficile de prescrire si le sol répondra en conditions drainées ou non-drainées. Pour ces situations, il est conseillé de mener deux calculs séparés, à la fois sous l’hypothèse de comportement drainé et non-drainé, et de définir la réponse globale via l’enveloppe de ces deux calculs.


Amortissement
 - Quand il s’agit de modéliser le comportement des couches de sol via une loi élastoplastique, l’amortissement matériel anticipé est obtenu directement via la dissipation qui résulte de la loi de comportement. Or, cette définition de l’amortissement diffère de la caractérisation conventionnelle de l’amortissement (accompagnant habituellement les définitions des courbes G-γ) via une courbe β-γ, où β représente le taux d’amortissement visqueux équivalent. Cette différence vient du fait que le type d’amortissement n’est pas le même dans les deux cas : d’une côté il y a l’amortissement visqueux β, d’autre part il y a l’amortissement hystérétique de la loi élastoplastique. Il s’ensuit de cette différence de définition qu’il n’est pas commode de faire un calage du modèle non-linéaire sur la base de la courbe β~γ. Il est toutefois important de vérifier a posteriori que les taux d’amortissement engendrés par le calage du modèle élastoplastique ne sont pas très différents par rapport aux valeurs de β qui correspondent aux courbes β-γ théoriques. En dehors de l’amortissement hystérétique, issu directement de la loi élastoplastique, il est également nécessaire de définir un amortissement « élastique », c’est-à-dire un taux d’amortissement de base, à appliquer dès le début de la sollicitation. La définition de cette composante d’amortissement (qui est calée à 1% - 3%) sert pour le bon conditionnement de l’équation d’équilibre dynamique et peut être introduite via un amortissement de type Rayleigh. Enfin, il est possible d’introduire une troisième composante d’amortissement, celle de l’amortissement numérique, qui résulte du choix d’un schéma d’intégration dissipatif. Cette troisième composante doit être utilisée avec précaution ; pour des applications courantes, elle est à éviter.


Liquéfaction
 - Un calage tout particulier des modèles de sol est exigé dans le contexte des études du risque de liquéfaction. Dans ce cas, il est nécessaire de faire appel à un modèle bi-phasique pour les zones potentiellement liquéfiables. Ceci permet de modéliser la montée de pression interstitielle qui survient dans ces matériaux peu perméables, lors d’une secousse sismique, phénomène qui peut conduire à l’état de liquéfaction. Le calage effectué dans ce cas consiste à montrer qu’un état liquéfié est bien obtenu sur une éprouvette de sol, consolidée à l’état de contrainte de la zone potentiellement liquéfiable, lorsqu’elle est sollicitée par le séisme théorique conduisant à la liquéfaction. Les paramètres les plus pertinents pour le calage vis-à-vis de la liquéfaction sont ceux qui gouvernent la réponse du sol en chargement selon l’axe hydrostatique (paramètres volumétriques), comme par exemple, l’angle de dilatance ou la règle de définition des déformations volumétriques plastiques etc.

 
5) Conditions aux limites : Frontières absorbantes - Modèles flottants


Dimensions de la zone à mailler
 - Un autre point subtil concernant la mise en œuvre d’une approche directe, concerne le « dimensionnement » du maillage et la définition des conditions aux limites. Tout d’abord, les dimensions du maillage doivent être suffisantes pour la description des phénomènes en jeu. Dans la direction verticale, le maillage est habituellement étendu jusqu’au niveau du substratum, où il peut être supposé que le comportement reste linéaire. Dans la direction horizontale, le maillage est étendu à une distance caractéristique égale à 4 à 5B, de part et d’autre de la structure, où B est la dimension caractéristique de la structure. Il s’ensuit que la dimension minimale d’un maillage pour la mise en œuvre d’une méthode directe d’IDSS est 10 fois la dimension caractéristique de l’ouvrage. Il est compris que ce critère, combiné avec les critères pour la taille des mailles, présentés plus haut, peut rendre prohibitif le coût d’un calcul d’IDSS par la méthode directe.


Conditions aux limites
 - Pour les problèmes où la sollicitation extérieure vient de la structure, la définition des conditions aux limites consiste à introduire des éléments de frontière absorbante le long du contour du maillage, dont le rôle est d’annihiler les réflexions parasites des ondes qui arrivent jusqu’aux extrémités du maillage. Les éléments de frontière absorbante couramment utilisés dans la pratique sont des amortisseurs visqueux standards qui ne sont exacts (en tant que frontières absorbantes) que pour une direction de propagation et une fréquence d’onde données. Pour passer outre ces limitations, plusieurs formulations de frontières absorbantes améliorées sont disponibles. Néanmoins, leur mise en œuvre reste relativement lourde et ces formulations restent peu utilisées dans les études courantes.

Pour les problèmes sismiques, le rôle des frontières absorbantes est à la fois la réduction des réflexions parasites et aussi l’introduction du mouvement sismique incident le long du contour du domaine modélisé en éléments finis. Pour atteindre ce deuxième objectif, il est souvent nécessaire de déterminer le mouvement au niveau du substratum (car en règle générale, les normes parasismiques fournissent le séisme de dimensionnement en surface du sol) : ce mouvement est utilisé comme le mouvement incident à la base du modèle. De plus, il faut déterminer le mouvement (en conditions de champ libre) sur une colonne de sol (du substratum jusque la surface), ce qui servira pour établir les conditions aux limites le long des parois latérales du modèle. Par ailleurs, la nécessité d’avoir un maillage suffisamment étendu dans la direction horizontale est expliquée par le fait que nous souhaitons établir le mouvement du champ libre à ses extrémités latérales. Pour les maillages dissymétriques, le calcul de la colonne du sol en champ libre doit être répété séparément pour la colonne du côté gauche et du côté droit du maillage.

Un contrôle qui doit toujours être fait, consiste à vérifier que le mouvement calculé en un nœud proche des frontières latérales du modèle est identique ou très proche du mouvement en champ libre. Cette vérification permet de s’assurer que le maillage est bien dimensionné pour le problème d’interaction étudié.


Phasage
 - Comme pour l’ensemble de problèmes géotechniques, il est aussi nécessaire pour les problèmes dynamiques, d’établir le vrai état de contrainte dans les géomatériaux avant de pouvoir introduire le chargement dynamique. Habituellement, toute étude d’IDSS est initiée par une phase de mise en gravité. A la fin de cette phase, les déplacements du système sont mis à zéro et les états de contrainte sont retenus dans la totalité du maillage. Suite à la complétion de la mise en gravité, il est possible d’introduire le chargement dynamique dans le système.


Modèles flottants
 - Toujours pour les problèmes sismiques, il est aussi important que la modélisation adoptée permette d’injecter simultanément la composante horizontale et la composante verticale du mouvement sismique. Pour la prise en compte de la composante verticale, on est souvent obligé à recourir à la méthode, dite « de modèle flottant » qui consiste à effectuer les étapes suivantes :

 
6) Choix d’un schéma d’intégration

La mise en œuvre de la méthode directe pour l’IDSS est complétée avec la définition du schéma d’intégration temporelle de l’équation d’équilibre dynamique. Les principes qui s’appliquent aux problèmes dynamiques non-linéaires généraux sont également applicables dans le cas des problèmes d’IDSS. Comme, il s’agit de problèmes non-linéaires, le schéma d’intégration temporelle et couplé avec une méthode d’intégration des lois de comportement non-linéaires, typiquement la méthode de Newton-Raphson ou la méthode de Newton-Raphson modifiée etc.

En règle générale, les problèmes de très courte durée (chocs, explosions etc.) sont traités par des schémas explicites. Ces schémas sont très efficaces si la matrice de masse peut être triangulée mais ils sont conditionnellement stables (les critères de stabilité conduisent souvent à des pas de temps très faibles). Pour les problèmes sismiques, on fait appel à des schémas implicites, le plus couramment utilisé étant le schéma d’intégration de Newmark (qui est inconditionnellement stable en régime linéaire pour un jeu donné des paramètres d’intégration).

Le choix d’un schéma dissipatif (entraînant de l’amortissement numérique) peut également être proposé de manière à faciliter la convergence de la solution pour des problèmes fortement non-linéaires.

Dans tous les cas, le choix d’un schéma d’intégration, approprié à un problème donné, avec la définition de tous les paramètres associés (tolérances, nombre d’itérations, stratégie pour gérer la divergence etc.) nécessite une très bonne connaissance des notions théoriques et aussi une expérience importante quant à la typologie des problèmes traités (sismiques, chocs, dynamique rapide, problèmes vibratoires etc.). Comme principes généraux, il est conseillé de commencer les analyses avec des schémas non-dissipatifs et avec des tolérances relativement strictes (de l’ordre de 10-4). A la fin de chaque analyse, il faut examiner avec prudence le fichier des résultats détaillés afin de vérifier que les résidus de chaque incrément ne dépassent pas les tolérances et que la convergence est obtenue sur la totalité des mailles du modèle.