ntroduction

La littérature mondiale a produit une quantité phénoménale d'articles traitant des éléments finis. De nombreux livres expliquent cette méthode née expérimentalement chez le constructeur d’avions Loockeed dans les années 1950. C’est en 1968 que les fondements mathématiques des éléments finis ont été posés et ils sont maintenant utilisés dans de nombreux programmes de calcul, l’évolution des moyens informatiques facilitant leur usage.

Les débuts de l'usage de cette méthode sont analysés par Radhakrishnan et al (1970) et depuis, les différents articles ont portés sur des développements de lois de comportement de plus en plus sophistiquées, ainsi que sur les techniques numériques de mise en œuvre (méthodes de résolution, mailleurs et post-processeurs).

En 1971, Duncan trouve que le développement des éléments finis n'est pas assez rapide et ceci essentiellement pour des raisons de prix. Il estime qu'il faut quelques mois à un ingénieur pour modifier correctement un programme (quelques semaines pour s'en servir) et que les crédits machine ne sont jamais suffisants. Si les crédits machine ne sont plus une limite, la machine peut en être une, la maîtrise de l’usage de ces outils nécessite encore beaucoup de temps.

Cette annexe a pour but de donner quelques repères pour la compréhension de cette méthode complexe.

Introduction aux éléments finis

3 - 1 - Formulation vectorielle

Le problème physique peut s'énoncer :

- Soit un solide V fixé sur sa frontière, dans le champ de pesanteur, il est soumis à des forces de volume et des forces extérieures en certains points. Il en résulte un déplacement .

Du point de vue formalisme mathématique et sont des vecteurs de dimension infinie (il y a une infinité de points dans le domaine ) et si est connu, est l'inconnue qui appartient à l'espace , espace des déplacements du solide V.

La transformation : ,est une application de l'espace vers un espace .

Cette transformation permet d'obtenir facilement l'énergie élastique (à travers la loi de comportement) du solide V que l'on notera Le travail des forces extérieures qui est défini par un produit scalaire est noté Le principe de l'énergie minimale nous donne comme solution de :

w est une forme bilinéaire et .est linéaire en., d'où l'écriture équivalente :

.

Ce problème de minimum est insoluble car on ne sait pas travailler sur des espaces infinis. Pour le rendre fini, on choisit est de dimension 2n et il existe une base 2n avec :

Le problème équivalent devient :

en posant

,

il s'agit de minimiser ce qui revient à résoudre :

système qui possède une solution car K, matrice de rigidité, est définie positive, coercive et symétrique.

La résolution du système sera facilitée si K est une matrice la plus creuse possible.

L'idéal serait d'avoir En fait on obtient une matrice bande, dont le choix des est difficile. On choisit pour les des supports dont l'intersection est la plus petite possible en regroupant quelques éléments d'un pavage d'éléments finis (Armand, Frémond, 1975).

3 - 2 - Approche classique

3 - 2 - 1 - Formulation élémentaire

Le milieu continu est idéalisé en un ensemble d'éléments (triangles dans cet exemple) dont on va écrire l'équilibre. L'énergie extérieure est produite par les déplacements des sommets, ce qui revient à dire que les forces sont concentrées à ces sommets appelés aussi noeuds.

A l'intérieur de ces éléments, une loi simple régit les déplacements. Comme on ne veut connaître les variables qu'aux noeuds, on va donner l'expression de l'énergie élastique et du travail des forces, en fonction du déplacement des noeuds.

Pour un triangle à trois noeuds, on note dans le cas plan :

les forces appliquées aux noeuds 1, 2, 3, et

les déplacements de ces mêmes noeuds.

Le travail des forces extérieures est alors

Dans le triangle, on exprime la loi de déplacement en chaque point de coordonnées x, y sous la forme linéaire suivante :

.

Les déformations sont alors :

, , .

Si on pose : on a formellement :

ou encore :

Le solide étant supposé élastique dans cette étape, les équations de Lamé s’écrivant :

l'énergie élastique de déformation s'écrit :

,

représente l'élément et est appelé aussi énergie interne.

Cette intégration est réalisable numériquement. Dans le cas du triangle à trois nœuds, elle est triviale car est une constante. Pour l'élément d'épaisseur h et de surface S on a :

L'énergie élastique est une intégrale de volume, mais ici h peut être pris égal à l'unité.

Les énergies mises en œuvre (extérieure et intérieure élastique) sont égales et ce pour tout déplacement licite, donc :

ou :

soit :

avec :

Par la matrice , que l'on découvre symétrique, on a établi une relation entre les déplacements aux nœuds et les forces aux nœuds, et ce pour un élément.

3 - 2 - 2 - Assemblage des éléments

Réunir les éléments dans la structure, c'est écrire une seule relation entre les forces appliquées à tous les n nœuds :

.et les déplacements à ces n nœuds :

Pour l'élément k de sommets l, m, p, dans la discrétisation du milieu continu, on a :

avec :

lème colonne pième colonne

mième colonne

donc :

Le déplacement est continu sur tout le milieu, vu la continuité aux nœuds. L'énergie externe s'écrit :

et l'énergie interne :

et ceci pour tout , donc :

avec

est la matrice de rigidité globale dont la construction nécessite une somme de produits de matrices.

Résolution du système obtenu

Le plus fréquemment c'est la méthode de Gauss qui est utilisée, avec des améliorations qui tiennent compte de l'allure particulière de la matrice obtenue.

On a ainsi les méthodes : Gauss symétrique, Cholewski et Gauss bande.

Si nous appelons n la dimension de la matrice carrée et m la demi largeur de la bande et que nous fassions le décompte des opérations, nous obtenons le tableau suivant :

Résolution du système obtenu

Le plus fréquemment c'est la méthode de Gauss qui est utilisée, avec des améliorations qui tiennent compte de l'allure particulière de la matrice obtenue.

On a ainsi les méthodes : Gauss symétrique, Cholewski et Gauss bande.

Si nous appelons n la dimension de la matrice carrée et m la demi largeur de la bande et que nous fassions le décompte des opérations, nous obtenons le tableau suivant :

Opérations

     

+

     

     

:

     

     

Nbre total

d'opérations

Gauss

cas général :

     

     

     

          

Cholewski

     

     

     

     

n

     

Gauss

symétrique

     

     

     

          

Gauss

matrice bande

     

          

          

La méthode utilisée donne un aperçu de la rapidité d'un programme. Il y a encore années le nombre d'opérations n'était pas le paramètre clé ; il fallait aussi envisager l'occupation de la mémoire de l'ordinateur. Avec les mémoires vives de grande taille et les gestionnaires de mémoire virtuelle, l’utilisateur n’a plus à se préoccuper de cette gestion. On trouvera dans Faure (1982) des éléments pour utiliser des mémoires réduites. Le paragraphe ci-après résume ces travaux pour utiliser de grandes matrices dans peu de mémoire vive.

3 - 4 - La matrice de rigidité globale et son implantation en machine.

La matrice dont nous avons décrit la constitution, est une matrice bande, symétrique. Elle peut donc être stockée sous la forme d'un tableau n X m , m étant la demi largeur de bande et n la dimension du système.

Même sous cette forme on atteint très vite les capacités maximales des ordinateurs des années 1980.

Divers algorithmes ont été proposés (par exemple Akras et Dhatt (1976)) pour minimiser cette largeur de bande.

Un stockage de la bande peut être fait en considérant son aspect sous forme d'un histogramme - Sky line Storage - (Thompson et al. 1980), ce qui permet dans certains cas de gagner encore un peu de place ; on ne se limite plus à la largeur maximale de bande. Mais l'usage de la mémoire annexe des ordinateurs donne d’autres possibilités.

Comme la mémoire centrale, celle où s'effectuent les calculs, était chère et limitée, on ne plaçait dans cette mémoire que ce qui était nécessaire aux calculs et le reste des informations était stocké sur des mémoires périphériques (disques).

Il fallait alors rajouter, à l'algorithme de calcul, un algorithme de gestion pour commander les transferts entre la mémoire centrale et la mémoire annexe.

Recuero et al (1979) découpent ainsi une méthode de Gauss jusqu'à ne se servir que de deux lignes du système à résoudre. On pouvait alors se servir des micro-ordinateurs qui arrivaient sur le marché.

Strabowski (1981) décrit un usage très élaboré de ce type d'utilisation. Fengels et Troost (1980) cherchent à diminuer le temps de transfert en limitant le nombre de transferts. Ils utilisent alors une partie de la mémoire centrale comme zone de stockage et ce n'est que lorsque cette zone est saturée qu'ils effectuent le transfert. C'est la notion de "buffer" des informaticiens qui est utilisée dans cet algorithme.

En poussant à l'extrême ces considérations de découpage, Verruitj (1980) propose un calcul d'éléments finis pour 100 éléments sur un micro calculateur (PET de Commodore), avec seulement 8000 octets de mémoire centrale !

3 - 5 - La méthode frontale

Cette méthode évite la construction de la matrice de rigidité globale en avançant pas à pas (élément par élément) dans la structure et en écrivant les relations avec seules les variables utiles.

Avec une telle approche, l'écriture, avec usage des mémoires auxiliaires, devient naturelle. C’est une solution élégante au problème de faible mémoire. (Faure, 1982). Aujourd’hui ce sont des blocs, ensembles d’éléments qui sont traités de cette façon avec parfois pour chaque blocs des résolutions faisant appel à différentes méthodes.

     

          

La non-linéarité

L'élasticité n'est pas le comportement réel du sol qui est bien plus compliqué. L'idée naturelle est de décomposer ce comportement quelconque en parties (ou domaines) où l'algèbre linéaire peut être utilisée. Ce sont alors des applications tangentes qui vont être décrites, la difficulté étant d'assurer une continuité entre ces domaines décrits linéairement.

En fait, la matrice d'élasticité d'un domaine est aussi fonction de l'histoire du matériau dans ses étapes précédentes et il est nécessaire d'élaborer des modes de mémorisation.

Mais comment déduire de tels modes à partir de l'expérimentation ?

Un article d'Argiris et al (1979) répertorie un certain nombre de méthodes non linéaires.

4 - 1 - Aspects de la non linéarité

4 - 1 - 1 - La notion de critère

La notion de critère a été introduit en plasticité, c'est une limitation de la valeur des contraintes. Pour que cette limitation soit facilement maniable, elle est donnée sous forme d'une fonction dont le signe indique la non plasticité du milieu. Pour une valeur nulle du critère, nous avons le seuil avec le début du comportement plastique.

Un très grand nombre de critères ont été proposés. Ces critères sont des descriptions de surface dans un espace de contraintes.

A la notion de critère, Hill ajoute la notion de normalité, (matériau standard) issue du principe de travail maximal, qui permet de différencier la déformation plastique de celle qui est élastique. La formulation s'alourdit, mais permet le calcul de la surcontrainte que l'on redistribue. On peut ainsi itérer sans modifier la matrice de rigidité, mais pas sans travail.

Explicitons un peu cette forme de calcul.

Lorsque le critère de plasticité est satisfait, le point représentatif atteint la surface de l'espace qui représente le critère, il se produit un écoulement plastique. L'accroissement total de déformation comprend deux parties, l'une élastique, l'autre plastique et l'on écrit :

.

Le principe du travail maximal établi par Hill permet d'affirmer que la surface représentant le critère est convexe et que est dirigé suivant une normale extérieure à cette surface

On en déduit que :

Dans le cas des sols, en tenant compte de l'écrouissage, f s'écrit et f varie au fur et à mesure des déformations plastiques (l'écrouissage est représenté par k).

Posons ces relations sous forme matricielle, facilement représentables dans un ordinateur.

.

A l'équilibre (sur la surface représentant le critère) on a ou en différenciant :

En petite déformation l'usage est de remplacer par e, ce qui s'écrit :

étant la matrice de rigidité

donc

en posant : où A apparaît comme un paramètre lié à l'écrouissage, on a :

on a et soit

et en reportant dans l'équation précédemment écrite, on obtient :

ou

On passe de appliquée dans l'élément, aux forces correspondantes (forces de reéquilibrage) appliquées aux nœuds de l'élément et l'on résout :

Exemple

Pour illustrer cette manipulation de formule, plaçons-nous en déformation plane dans le cas du critère de Coulomb. Ce dernier s'écrit en fonction de deux contraintes principales seulement : :

On identifie A et k en écrivant le critère sous la forme :

avec

il vient :

avec

4 - 3 - Problèmes posés par les méthodes de résolution

La non linéarité se traduit par la définition d'un critère. C'est une surface dans l'espace des contraintes et une loi d'écoulement, dépendante du critère, qui est décrite dans les programmes. Cette loi fournit alors un à redistribuer et les ajustements successifs sont arrêtés quand la surcharge devient inférieure à un seuil fixé.

Les différentes méthodes de redistribution ne sont pas très convaincantes et comportent le risque de mauvaise convergence. Cette convergence serait plus efficace dans un schéma Newton, où l'on recalcule à chaque fois la matrice tangente, alors que l'usage veut l'emploi d'une méthode Newton modifiée où la matrice tangente est celle de départ.

Méthode de NEWTON

     

     

Méthode de NEWTON modifiée

Une "bonne" approche serait de se débarrasser de ces redistributions et donc de ne pas trop s'éloigner de la courbe effort déformation, c'est-à-dire travailler en matrice tangente, avec des pas de chargement petits, mais que veut dire ‘petit’.

Pour illustrer ces difficultés de convergence examinons deux méthodes d’itérations..

4 - 3 - 1 - La méthode des sous-incréments.

Pour présenter la méthode des sous-incréments, Bushnell (1977) commence par rappeler la méthode usuelle et fait la différence entre deux niveaux d'itération (ou boucles) notés A et B.

    * Méthode usuelle

Pour un incrément de charge on résout un système A de n équations algébriques dont le rang est égal au nombre de degré de liberté du système modélisé. On détermine alors la déformation plastique et la direction de l'écoulement plastique est supposée constante dans tout l'incrément (loi de normalité en général). Ceci est résolu par itération dans un système "local" d'équation, au niveau de l'élément, et on appelle B ce système.

Tillerson montre que la convergence n'est pas toujours assurée du fait que les charges et décharges entre chaque itération, sont appliquées à travers une matrice A qui varie de façon discontinue.

    * Méthode des sous-incréments

Elle se différencie de la méthode usuelle par les item suivants :

    - le calcul de ne se fait pas dans la boucle A et on peut faire varier la direction de l'écoulement.

    - la convergence de la redistribution des "charges" plastiques est assurée.

Dans le processus des sous-incréments ; l'incrément total de déformation entre deux cas de charge est divisé en sous-incréments.

Pour chaque sous-incrément, la direction de l'écoulement plastique est prise constante et est normale à la surface du critère, en un point correspondant au précédent sous- incrément. Et pour chaque sous-incrément de déformation, le sous-incrément correspondant de contrainte est déterminé par la loi d'écoulement et la loi effort- déformation (uniaxial stress strain curve).

Ainsi le système B est résolu pour chaque sous-incrément et en chaque point en s'appuyant directement sur une loi expérimentale.

Cela revient à dire que dans l'incrément de charge, on raisonne avec les déformations que l'on s'impose pas à pas.

4 - 3 - 2 - La méthode ILI (Integrated Load Increment)

Ce processus itératif est présenté par Cheng et al (1980) comme une extension de la méthode de (Zienkiewicz et al, 1969). Une économie de calcul est faite en intégrant préalablement la loi effort-déformation sur l'intervalle (pas de chargement), ce qui évite la redistribution des efforts dans le cycle suivant.

Les séquences sont dans l'ordre.

a/ Calcul de la matrice d'élasticité élémentaire de l'élément :

.B matrice de dérivation donnant les déformations à partir des déplacements.

. matrice d'élasticité.

b/ Application de l'incrément de charge

c/ Calcul de .1ère itération

d/ Calcul de matrice de plasticité.

e/ Calcul de (1ère itération).

f/ Calcul de .

g/ Calcul de par la formule obtenue par intégration de la loi d'écoulement.

h) Calcul de forces non compensées par :

où K est la matrice tangente.

i/ Refaire les étapes c à h jusqu'à ce que soit inférieur à une valeur définie à l'avance.

j/ Calcul de toutes les valeurs de fin d'étapes.

k/ Passage au chargement suivant.

La figure suivante est une représentation schématique de la méthode ILI.

Représentation schématique de la méthode ILI

(d'après Cheng et Al. 1980)

On peut économiser l'étape c en extrapolant à partir de l'itération précédente.

Le gain en temps calcul est à l'étape g, car il n'y a pas d'itération pour obtenir .

4 - 3 - 3 - Résolution par accroissement des déplacements

Quand on arrive près de la charge limite, il peut être difficile, voire impossible, de redistribuer les forces non compensées si l'incrément de charge Q est trop important ; ce qui se traduit dans certains algorithmes par la prise en compte de Q décroissants.

Une autre idée est alors de travailler par accroissement de déplacement, puis de calculer la charge correspondante.

Cette idée apparaît en 1977 (Walter et Al.). Battoz et Dhatt (1979) décrivent un algorithme que Pauwell et Simon (1981) raffinent. Le champ de recherche de ces auteurs est le flambement des structures et ils trouvent ainsi un moyen d'éviter la notion de bifurcation par un suivi très réaliste et facile d'une courbe avec asymptote horizontale. Cependant, ils traitent le déplacement point par point, ce qui est admissible dans le cas d'une structure élancée qui flambe ; mais paraît assez délicat pour être mis en œuvre pour les sols : quel déplacement en quel point ?

Cette voie n'est cependant pas à rejeter, elle pourrait être employée pour les problèmes d'ancrage où l'on mesure le déplacement, mais pas la contrainte.

5 - Conclusion

Traiter des comportements non linéaires en éléments finis n’est pas une chose numériquement et algorithmiquement simple. Si les programmes des années 90 comportent des algorithmes encore plus élaborés que ceux présentés ci-dessus, malgré les efforts des auteurs, leur robustesse n’est pas absolue et des cas particuliers de chargement peuvent les mettre en défaut. L’art de l’ingénieur est de détecter et contourner ces résultats erronés

Loi de Duncan, programme FEADAM.

Développé pour le calcul des barrages en terre construits par remblaiement hydraulique, le modèle hyperbolique de Duncan correspond aux sols contractants. L’idée est de suivre pas à pas le chargement, découpé en incréments, en évitant les décharges importantes.

Pour un pas de chargement, la loi de Hooke, l’élasticité, s’écrit :

Pour obtenir les deux paramètres de l’élasticité, ici E et B, on se sert de la ressemblance des courbes effort-déformation d’un sable lâche avec une hyperbole. On écrit :

(s 1 - s 3) = e / ((1 / Ei) + e / (s 1 - s 3)ult))

Ei est le module tangent initial, (s 1 - s 3)ult la valeur asymptotique qui permet à l’hyperbole de ressembler le mieux possible à la courbe expérimentale. Cette dernière présente un palier que l’on note (s 1 - s 3)rupt. Et l’on note Rf = (s 1 - s 3)rupt / (s 1 - s 3)ult qui apparaît comme une constante pour un sol donné. On a 0.5 < Rf < 0.9

Janbu, suite à de nombreux essai, détermine Ei suivant :

Ei = K pa (s 3 / pa)n

avec K,n, paramètre et exposant du module, pa pression atmosphérique pour rendre l’expression sans dimension.

La loi de Coulomb s’écrit :

(s 1 - s 3)rupt = (2c cos j + 2 s 3 sin j ) / (1 - sin j )

On a donc relié s 1, s 3, avec e . Ce que l’on cherche c’est Et, le module tangent. Par dérivation on obtient :

Et = ( 1 - (Rf (1 - sin j ) (s 1 - s 3) ) / (2c cos j + 2 s 3 sin j ) )2 K pa (s 3 / pa)n

Ceci permet de simuler un chargement continu. Pour un déchargement, en négligeant l’hystérésis un seul module est nécessaire et par analogie on écrit :

Edr = Kdr pa (s 3 / pa)n

L’exposant n est pris égal à celui en chargement et Kdr apparaît supérieur à K. (de 20% pour un sable serré à 300% pour un sable très lâche. )

B, le second terme du modèle, est indépendant de (s 1 - s 3). Il est défini par :

B = Kb pa (s 3 / pa)m

avec 0 < m < 1 et Et / 3 <B < 17 Et. En effet, si B = Et / 3, la valeur du coefficient de Poisson,

n t = 0.5 - Et / 6 B

devient nulle, ce qui est à éviter numériquement et physiquement. De même pour ne pas avoir n > 0.5 on limite B à 17 Et.

La loi de Duncan intègre aussi le fait que l’angle j dépend de la pression de confinement s 3, lorsque celle-ci varie de façon importante, ce qui est le cas pour un barrage. Il écrit :

j = j 0 - D j log ( s 3 / pa)

La loi de Duncan est donc définie par neuf paramètres : c, j 0, D j , K, Kb , Kdr , n ,m, Rf.

Pour déterminer ces neufs paramètres toute une procédure, à partir d’essais triaxiaux a été proposée. (Duncan et Chang, 19XX). Les résultats sur de nombreux sols ont aussi été tabulés.

Modèle de Duncan

L'élasticité.

Au début, il y avait l’élasticité, voyons comment ceci a conduit aux éléments finis.

L’élasticité a été énoncée par Hooke en 1678 sous la forme de l’anagramme :

c e i i i n o s s s t t u v explicité deux ans après par la formule latine : ut tensio sic vis

Il y a donc proportionnalité entre l'effort et la déformation, ce qui montre une relation biunivoque entre ces deux grandeurs. Ceci implique que sans effort la déformation est nulle, il y a retour à l'état initial. La relation est linéaire et l'élasticité est linéaire C'est celle que nos ordinateurs acceptent le plus facilement, l'algèbre linéaire étant la base de nos calculs.

La formule la plus générale de l'élasticité linéaire reliant les tenseurs de contraintes et de déformations :

Avec la symétrie des tenseurs , et avec les considérations thermodynamiques (énergie dissipée nulle), les 81 composants du tenseur se réduisent à 21 valeurs indépendantes. Les propriétés de symétrie spatiale de la matière permettent de réduire encore ce nombre.

2 - 1 - Différentes élasticités

    * Elasticité orthotrope.

    Elle se caractérise par une symétrie par rapport à trois plans orthogonaux pris 2 à 2.

    On obtient dans ce cas la relation :

    

ce qui fait 9 paramètres pour définir l'orthotropie.

La détermination expérimentale de G, H, J (essai de cisaillement homogène) est bien plus délicate que celle de A, B, C, D, E, F. Ceci annonce les difficultés expérimentales d’imaginer des essais capables de fournir les paramètres de lois plus complexes.

    * Elasticité orthotrope de révolution ou aléotropie

Un sol est généralement stratifié horizontalement du fait de son histoire géologique, la direction verticale est privilégiée alors que deux directions horizontales perpendiculaires ne sont pas différentiables. On en déduit :

5 paramètres suffisent pour déterminer cette loi de comportement élastique.

    * Le modèle isotrope

N'importe quel axe de l'espace est de révolution. On obtient alors :

et si l'on pose : appelés coefficients de LAME, on a :

Les coefficients de Lamé sont liés à E et n par :

Une étude de Cordary et al (1976) montre que les insuffisances du modèle isotrope défini par seulement deux paramètres peuvent être expliquées en partie par un modèle aléotrope défini par cinq paramètres. Car la difficulté est alors dans l'obtention de ces cinq paramètres. En détaillant l'approche, il apparaît que la détermination des paramètres isotropes est ambiguë car l'expérimentation se fait sur un échantillon de sol, aélotrope par nature, et les paramètres mesurés ne devraient donc pas être utilisés dans des calculs bâtis sur une loi de comportement isotrope.

Cette double erreur (utiliser des paramètres d’une loi dans une autre et mal définir les paramètres) peut expliquer à la fois les insuffisances de l'isotropie appliquée au sol et aussi les bons résultats obtenus parfois, les deux erreurs se compensant dans ce cas.

Les moyens de calcul se prêtent particulièrement bien aux formes linéaires. L'élasticité linéaire décrite ci-dessus n'est pratiquement jamais vérifiée dans la nature (il y a dissipation d'énergie) et les lois de comportement des sols, qu'il faut décrire sont éminemment non linéaires. Pour tourner cette difficulté et se ramener à une description linéaire, des lois tangentes, valables pour un certain incrément, vont être utilisées.

2 - 2 - Formulation du travail interne dans un solide élastique.

Soit un solide de frontière sur laquelle nous avons des actions de contact ,et soit un champ de déplacement licite.

Le travail de ces actions dans le champ de déplacement est :

en développant le tenseur des contraintes suivant sa définition : et à l'aide de la formule de Green on transforme l'intégrale de surface en intégrale de volume :

On sait que le tenseur des contraintes est symétrique (Cauchy) :

On introduit l'équation indéfinie du mouvement : pour avoir

ou

ce qui montre que l'énergie cinétique du système est égale au travail des forces extérieures diminué du travail des forces intérieures dont l'expression est le produit des tenseurs de contraintes et de déformations.