Web Analytics Made Easy - Statcounter
Privacy Policy Cookie Policy Terms and Conditions

[HOME PAGE] [STORES] [CLASSICISTRANIERI.COM] [FOTO] [YOUTUBE CHANNEL]


Méthode des éléments finis

Méthode des éléments finis

Page d'aide sur l'homonymie Pour les articles homonymes, voir Élément.
Solution bidimensionnelle d'une équation magnétostatique obtenue par éléments finis (les lignes donnent la direction du champ et la couleur son intensité)
Maillage utilisé pour l'image du haut (le maillage est plus resserré autour de la zone d'intérêt)
Simulation numérique d'un essai de choc sur une voiture  : les cellules utilisées pour le maillage sont visibles sur la surface du véhicule.

En analyse numérique, la méthode des éléments finis (MÉF, ou FEM pour finite elements method) est utilisée pour résoudre numériquement des équations aux dérivées partielles. Celles-ci peuvent par exemple représenter analytiquement le comportement dynamique de certains systèmes physiques (mécaniques, thermodynamiques, acoustiques, etc.).

Concrètement, cela permet par exemple de calculer numériquement le comportement d'objets même très complexes, à condition qu'ils soient continus et décrits par une équation aux dérivées partielles linéaire : mouvement d'une corde secouée par l'un de ses bouts, comportement d'un fluide arrivant à grande vitesse sur un obstacle, déformation d'une structure métallique, etc.

Introduction

La méthode des éléments finis fait partie des outils de mathématiques appliquées. Il s'agit de mettre en place, à l'aide des principes hérités de la formulation variationnelle ou formulation faible, un algorithme discret mathématique permettant de rechercher une solution approchée d’une équation aux dérivées partielles (ou ÉDP) sur un domaine compact avec conditions aux bords et/ou dans l'intérieur du compact. On parle couramment de conditions de type Dirichlet (valeurs aux bords) ou Neumann (gradients aux bords) ou de Robin (relation gradient/valeurs sur le bord).

Il s'agit donc avant tout de la résolution approchée d'un problème, où, grâce à la formulation variationnelle, les solutions du problème vérifient des conditions d'existence plus faibles que celles des solutions du problème de départ et où une discrétisation permet de trouver une solution approchée. Comme de nombreuses autres méthodes numériques, outre l'algorithme de résolution en soi, se posent les questions de qualité de la discrétisation :

  • existence de solutions,
  • unicité de la solution,
  • stabilité,
  • convergence,
  • et bien sûr : mesure d'erreur entre une solution discrète et une solution unique du problème initial.

La partie 2 va présenter le cadre général de la méthode des éléments finis, ainsi que le cas pratique le plus courant considérant des équations aux dérivées partielles linéaires dont on cherche une approximation par des fonctions affines.

La présentation en partie 3 est essentiellement physique, notamment mécanique. Elle ne doit être considérée que comme une présentation des éléments constitutifs de la modélisation discrète utilisée en résistance des matériaux via la méthode des éléments finis. C'est une approche tout à fait valide, un bon exemple pédagogique. Elle apporte un biais certain quant à une approche plus générale, du fait notamment de la linéarité supposée des matériaux.

Méthode des éléments finis

Principe général

l'utilisation de mailles tétraédriques (gauche) permet de mailler « fidèlement » des géométries complexes. L'utilisation d'un maillage régulier (droite) n'est possible que sur des domaines simples, mais permet de réduire le nombre de nœuds et donc le coût du calcul.

Considérons un domaine Ω (typiquement une portion de l'espace) dont la frontière est notée δΩ ou Σ. Nous cherchons à déterminer une fonction u définie sur Ω, qui est une solution d'une équation aux dérivées partielles (ÉDP) pour des conditions aux limites données. L'ÉDP décrit le comportement physique du système, il s'agit par exemple des lois de l'élasticité pour un problème de résistance des matériaux ou des équations de Maxwell pour les problèmes d'électromagnétisme. Les conditions aux limites sont les contraintes s'exerçant sur le système. Par exemple, pour un problème de résistance des matériaux, on impose le déplacement de certaines parties du système, par exemple, on impose qu'une zone d'appui soit immobile, et on impose des efforts sur d'autres zones (poids, pression de contact…).


La méthode des éléments finis (MÉF) permet de résoudre de manière discrète et approchée ce problème ; on cherche une solution approchée « suffisamment » fiable.

La discrétisation consiste à « découper » le domaine Ω, c'est-à-dire à chercher une solution du problème sur un domaine polygonal ou polyédrique par morceaux ; il y a donc une rédéfinition de la géométrie. Une fois la géométrie approchée, il faut choisir un espace d'approximation de la solution du problème, dans la MÉF, cet espace est défini à l'aide du maillage du domaine (ce qui explique aussi pourquoi il est nécessaire d'approcher la géométrie). Le maillage du domaine permet d'en définir un pavage dont les pavés sont les éléments finis.

Sur chacun des éléments finis, il est possible de linéariser l'ÉDP, c'est-à-dire de remplacer l'équation aux dérivées partielles par un système d'équations linéaires, par approximation. Ce système d'équations linéaires peut se décrire par une matrice ; il y a donc une matrice par élément fini. Cependant, les conditions aux frontières sont définies sur les frontières du système global et pas sur les frontières de chaque élément fini ; il est donc impossible de résoudre indépendamment chaque système. Les matrices sont donc réunies au sein d'une matrice globale. Le système d'équations linéaires global est résolu par l'ordinateur (des systèmes simples peuvent être résolus à la main et constituent en général des exercices d'apprentissage).

L'ÉDP est résolue aux nœuds du maillage, c'est-à-dire que la solution est calculée en des points donnés (résolution discrète) et non en chaque point du domaine Ω. Cela nécessite de pouvoir interpoler, c'est-à-dire déterminer les valeurs en tout point à partir des valeurs connues en certains points. On utilise en général des fonctions polynomiales.

Un élément fini est la donnée d'une cellule élémentaire et de fonctions de base de l'espace d'approximation dont le support est l'élément, et définies de manière à être interpolantes (voir Fonctions de base).

Nous voyons ici poindre trois sources d'erreur, c'est-à-dire d'écart entre la solution calculée et les valeurs réelles :

  • la modélisation de la réalité : le domaine Ω correspond en général à des pièces matérielles, le calcul se fonde sur des versions idéales (sans défaut) des pièces, de la matière et des conditions aux limites ; cette source d'erreur n'est pas spécifique à la méthode des éléments finis, et peut être prise en compte par la méthode contrainte-résistance ;
  • la géométrie idéale et continue est remplacée par une géométrie discrète, et les valeurs sont interpolées entre des points ; plus les points sont espacés, plus la fonction d'interpolation risque de s'écarter de la réalité, mais à l'inverse, un maillage trop fin conduit à des temps de calculs extrêmement longs et nécessite des ressources informatiques (en particulier mémoire vive) importante, il faut donc trouver un compromis entre coût du calcul et précision des résultats ;
  • s'agissant de calcul numériques, il se produit inévitablement des erreurs d'arrondi, les nombres étant représentés par un nombre fini d'octets.

Toute l'habileté de l'ingénieur consiste à maîtriser ces erreurs notamment :

  • en simplifiant la géométrie (defeaturing), en enlevant des détails qui se situent loin des zones que l'on veut étudier et ayant une faible influence sur le résultat ;
  • en choisissant des maillages adaptés, par exemple, des maillages de type poutre pour des pièces élancées, ou de type coque pour des pièces fines, en découpant la pièce pour pouvoir faire des maillages réguliers sur certaines zones, en affinant le maillage dans les zones critiques…
  • en ayant un regard critique sur le résultat.

Bien qu'il existe de nombreux logiciels exploitant cette méthode et permettant de « résoudre » des problèmes dans divers domaines, il est important que l'utilisateur ait une bonne idée de ce qu'il fait, notamment quant au choix du maillage et du type d'éléments qui doivent être adaptés au problème posé  : aucun logiciel ne fait tout pour l'utilisateur, et il faut toujours garder un œil critique vis-à-vis de solutions approchées. Pour cela il existe des indicateurs d'erreur et des estimateurs d'erreur qui permettent d'ajuster les différents paramètres.

La solution trouvée, il reste cependant à déterminer les caractéristiques de la méthode ainsi développée, notamment l'unicité de l'éventuelle solution ou encore la stabilité numérique du schéma de résolution. Il est essentiel de trouver une estimation juste de l'erreur liée à la discrétisation et montrer que la méthode ainsi écrite converge, c’est-à-dire que l'erreur tend vers 0 si la finesse du maillage tend elle aussi vers 0.

Dans le cas d'une ÉDP linéaire avec opérateur symétrique (comme l'est l'opérateur laplacien), il s'agit finalement de résoudre une équation algébrique linéaire, inversible dans le meilleur des cas.

Dimensions

Pour l'explication, nous développons ici la méthode des éléments finis en deux dimensions à valeurs réelles. On suppose que les équations étudiées sont des équations différentielles d'ordre deux.

La méthode est généralisable à des cadres d'espaces de dimension différente ou pour des équations aux dérivées partielles d'ordre supérieur :

  • on traite ici le cas d'une solution réelle à une ÉDP ; les cas où la dimension de la solution serait plus grande se traitent de façon similaire mais nécessitent des écritures plus complètes ; les cas les plus couramment rencontrés sont la dimension 1 (comme ici), 2 ou 3 (pour des problèmes de mécanique), 6 ou 12 (pour des problèmes d'électromagnétisme respectivement réels ou complexes) ;
  • les degrés de différentiation supérieurs sont ramenés à un degré moindre par la méthode classique de réduction de degré : on fait intervenir des variables supplémentaires, c'est-à-dire des dérivées partielles des variables de départ (exemple classique : les ÉDP de la mécanique statique des poutres font intervenir la dérivation partielle d'ordre 4) ; il est parfois possible, pour des degrés supérieurs, d'appliquer plusieurs fois les méthodes de formulation variationnelles afin d'obtenir des ordres plus faibles — en tout cas lorsque le degré de dérivation est pair.

Bien que théoriquement la méthode soit transposable en dimensions supérieures du support, techniquement la complexité de création des discrétisations croît avec la dimension… et pratiquement, on résout rarement des problèmes en dimensions supérieures à 3 — y compris des problèmes de dynamique en espace à 3 dimensions qui pourraient être traités en quatre dimensions mais sont traités en réalité avec une méthode mixte éléments finis « en espace » et en différences finies « en temps ».

Cadre algébrique, analytique et topologique

Soit un domaine (ouvert borné et connexe) Ω de \mathbb{R}^2, de bord δΩ, et d'adhérence (compacte) Ω. Pour simplifier les représentations, on suppose le bord polygonal.

Soient les fonctions de Ω dans \mathbb{R} différentiables sur Ω (compact) et deux fois différentiables sur Ω (ouvert). De telles fonctions sont continues et différentiables sur le bord du compact. Soit V(Ω) l'ensemble de ces fonctions (V est un espace vectoriel de dimension infinie et V0 est le sous-espace vectoriel de fonctions de V nulles sur le bord δΩ).

Soient les applications continues sur Ω et différentiables sur Ω, de carré sommables sur Ω et de gradient de carré sommable sur Ω (ou de dérivées partielles de carré sommable, ce qui revient au même avec le support de dimension finie). Nommons cet espace \mathcal{H}_1(\Omega). Cet espace est un espace de Sobolev. On dote cet espace vectoriel d'un produit scalaire issu de celui de L2 tel que si (u, v ) appartiennent à cet espace alors le produit scalaire de u et v est : \langle u | v \rangle_{ \mathcal{H}_1(\Omega) } = \int_\Omega{ ( \nabla u \cdot \nabla v + u v ) \,\mathrm{d}\omega }

On note \mathcal{H}_1^0(\Omega) le sous-espace vectoriel de \mathcal{H}_1(\Omega) dont les fonctions sont nulles sur le bord δΩ. L'opérateur (x,y) \mapsto x \cdot y est un produit scalaire sur l'espace \mathbb{R}^2

Cas organique

Hypothèses

On considère ƒ une fonction continue sur Ω de carré sommable et u la solution de l'équation aux dérivées partielles suivante sur Ω (\displaystyle \Delta est l'opérateur laplacien) : \displaystyle - \Delta u + k^2 u = f

Avec la condition au bord u = 0 sur δΩ. Ceci peut également se réécrire u ∈ V0. Cette condition au bord s'appelle la condition de Dirichlet.

On démontre qu'il existe une solution unique à ce problème d'EDP à l'aide du théorème de Lax-Milgram.

Formulation faible

Soit v ∈ V0 quelconque. Multiplions les deux parties de l'équation précédente par v puis sommons sur le domaine Ω, puisque v et ƒ sont tous deux de carré sommable sur ce domaine. On obtient l'équation : - \int_\Omega{ v \Delta u \,\mathrm{d}\omega } + k^2 \int_\Omega{ v u \,\mathrm{d}\omega } = \int_\Omega{v f \,\mathrm{d}\omega }\,

On utilise pour le premier terme une intégration par parties[1] : - \int_\Omega{ v \Delta u \,\mathrm{d}\omega } = - \int_{\partial \Omega }{ \frac{ \partial u }{ \partial n } v \,\mathrm{d}s } + \int_\Omega{ ( \nabla u \cdot \nabla v ) \,\mathrm{d}\omega }

Dans cette formulation, v est nulle sur le bord (v ∈ V0) ce qui permet d'obtenir la formulation faible du problème : \int_\Omega{ \nabla u \cdot \nabla v \,\mathrm{d}\omega } + k^2 \int_\Omega{ v u \,\mathrm{d}\omega } = \int_\Omega{ v f \,\mathrm{d}\omega }

Si u est deux fois différentiable, il y a équivalence entre cette formulation et celle du problème initial donné dans la section hypothèse et alors la solution de la formulation faible est la même que la solution initiale. On peut donc résoudre la formulation faible au lieu de résoudre le problème initial.

La question de savoir s'il y a équivalence entre la formulation faible et la formulation initiale donnée dans les hypothèses peut être particulièrement délicate dans les cas limites où l'ouvert Ω n'est pas suffisamment régulier (par exemple s'il y a des points singuliers) ou si ƒ n'est pas suffisamment dérivable (si l'on ne suppose pas que ƒ est au moins L^2). Il faut alors souvent se ramener à une étude au cas par cas et rien ne dit que la formulation faible aura les mêmes solutions que l'équation de départ. Dans la majorité des problèmes physiques la solution est souvent \mathcal{C}^\infty et l'on ne se pose pas de tels problèmes. Néanmoins, pour des domaines avec des points singuliers cette équivalence peut poser problème. Ceci peut être gênant pour l'étude de fissures en mécanique des milieux continus par exemple.

Notations et cadre général

Pour plus de généralité et pour rendre la suite plus lisible on utilisera les notations suivantes : a(u,v) = \int_\Omega{ \left(\nabla u \cdot \nabla v + k^2 u v \right) \,\mathrm{d}\omega } avec a un opérateur bilinéaire symétrique (de V2 dans \mathbb{R}) ; \mathcal{L}(v) = \int_\Omega{ f v \,\mathrm{d}\omega } avec \mathcal{L} un opérateur linéaire (de V dans \mathbb{R}).

On peut résoudre par la méthode des éléments finis toute équation aux dérivées partielles dont la forme faible se met sous la forme a( u, v ) = \mathcal{L}(v)

On montre que a est un opérateur bilinéaire coercif continu selon la norme \mathcal{H}_1^0 (cf. espace de Sobolev) et \mathcal{L} un opérateur linéaire continu également selon la norme \mathcal{H}_1^0. Avec ces notations, le problème se reformule ainsi : \forall v \in V_0, a( u, v ) = \mathcal{L}(v)

Puisque k>0, le théorème de Lax-Milgram assure l'existence de la solution et montre que u, solution du problème précédent, est la solution unique du problème d'optimisation de la fonctionnelle suivante : \forall v \in V_0, \mathcal{J}(v) = \frac{1}{2} a( v, v ) - \mathcal{L}(v)

Cette égalité peut avoir un sens physique notamment du point de vue de l'énergie pour certaines équations physiques et peut servir à montrer l'existence et l'unicité de la solution grâce aux propriétés de a et de \mathcal{L} (linéarité, coercivité, …).

Choix d'un maillage et discrétisation

Choix d'un maillage

Un exemple de maillage triangulaire
Raffinement de maillage par diminution de la taille d'éléments (type h, à gauche) ou par augmentation du degré des éléments (type p, à droite).

La méthode des éléments finis repose sur un découpage de l'espace selon un maillage. D'habitude l'on choisit un maillage carré ou triangulaire mais rien n'interdit de choisir des maillages plus complexes. Il n'est pas non plus nécessaire que le maillage soit régulier et l'on a tendance à resserrer le maillage près des endroits d'intérêt (par exemple aux endroits où l'on pense que la solution va beaucoup varier) ; cependant, il faut veiller à avoir des éléments faiblement distordus (se rapprocher d'un polygone régulier). Plus ce maillage est resserré, plus la solution que l'on obtient par la méthode des éléments finis sera précise et proche de la « vraie » solution de l'équation aux dérivés partielles.

On appelle traditionnellement h la plus grande dimension d'un élément (le diamètre de la sphère dans laquelle s'inscrit l'élément), et p le degré du polynôme décrivant le côté ou l'arête (p = 1 pour des côtés/arêtes droits, p = 2 pour des côtés/arête présentant une courbure).

Dans l'idéal, le maillage doit donc épouser les contours δΩ du domaine. Si δΩ est courbe, alors on peut :

  • soit utiliser des éléments plus petits, on parle de raffinement de type h ;
  • soit utiliser des éléments dont les côtés (en 2D) ou arêtes (en 3D) sont courbes, on parle de raffinement de type p.

Lorsque le degré des polynômes p est élevé, on parle de méthode des éléments spectraux.

Définition formelle

Le cadre mathématique des éléments finis permet de poser un système assurant la bonne définition de la solution.

Élément fini  On appelle élément fini la donnée d'un triplet (K, P_K, \Sigma_K) avec

  • K est un domaine géométrique,
  • P_K est un espace de fonctions sur K, qu'on appelle espace des fonctions de base,
  • \Sigma_K est un ensemble de formes linéaires sur P_K, qu'on appelle degrés de liberté.

Cette définition offre a priori beaucoup de liberté, mais en général, plusieurs conditions sont imposées : le domaine K est pris comme non dégénéré, l'espace des fonctions de base sera choisi de dimension finie et simples à calculer et les degrés de liberté seront pris comme vérifiant la propriété d'unisolvance :

Unisolvance  L'élément fini (K, P_K, \Sigma_K) est dit unisolvant si spécifier les valeurs sur chaque degré de liberté permet de spécifier une unique fonction de P_K.

Fonctions de base

Fonction de base en dimension 1. Les xi sont les nœuds du réseau.

On doit après prendre une base de fonctions « adaptées » au maillage. Plusieurs choix sont alors possibles. En général, les fonctions de base utilisées pour les éléments finis sont interpolantes, c'est-à-dire que les valeurs nodales sont les valeurs des grandeurs inconnues aux nœuds.

La plus simple est l'emploi des polynômes de Lagrange. Dans cette méthode les fonctions de base valent 1 à un nœud du maillage et 0 à tous les autres. La fonction de base i est alors la fonction valant 1 au nœud i et 0 sur les autres nœuds et qui est polynomiale sur chaque élément. Un exemple de telles fonctions est représenté en dimension 1 à côté. Il y a autant de fonctions de base par élément que de nombre de nœuds.

On appelle élément la donnée d'une géométrie (souvent polygonale en 2D, polyédrique en 3D) et de fonctions de base associées à cette géométrie.

D'autres solutions peuvent exister pour les fonctions de base. On cite ici un seul exemple les éléments finis d'Hermite qui ont la particularité d'avoir deux fonctions de base associées à chaque nœud. Dans cette version, la valeur de la solution est ajustée avec la première fonction alors que la deuxième permet d'ajuster la valeur de la dérivée. Ce type de fonctions de base peut avoir un intérêt pour la résolution de certaines équations aux dérivées partielles (par exemple l'équation des plaques en mécanique des milieux continus), même si elle nécessite d'avoir deux fois plus de fonctions pour un maillage donné.

Quelques éléments classiques
Principaux types d'éléments utilisé en 3D.

En 2D :

  • triangles
    • triangles de degré 1, (triangles à 3 nœuds, fonctions linéaires)
    • triangles de degré 2 (triangles à 6 nœuds, polynômes de degré 2)
  • quadrilatères
    • quadrilatères de degré 1 (carrés à quatre nœuds, fonctions linéaires)
    • quadrilatères de degré 2 (carrés à 8 ou 9 nœuds, polynômes de degré 2)

En 3D :

  • tétraèdres
    • tétraèdres de degré 1, (quatre nœuds, fonctions linéaires)[2]
    • tétraèdres de degré 2, (dix nœuds, polynômes de degré 2)[3]
  • hexaèdres[4]
    • hexaèdres de degré 1, (huit nœuds, fonctions linéaires)
    • hexaèdres de degré 2, (vingt nœuds, polynômes de degré 2)
    • hexaèdres triquadratique, (vingt-sept nœuds, polynômes de degré 2)


Les deux fonctions de base d'Hermite associées au nœud 0 avec les nœuds voisins en +1 et -1 en dimension 1
Discrétisation

Soit le maillage \mathcal{M} et la base b=(e_1...e_n) associée. Puisque la condition de Dirichlet impose des fonctions nulles aux bords, on utilise uniquement la sous-base b limitée aux points intérieurs de \Omega.

On cherche la solution \bar{u} du problème discrétisé ainsi :

 \bar{u} \in V_n^0 | \forall v \in V_n^0, a( \bar{u}, v ) = \mathcal{L}(v)

Or dans cet espace discrétisé, dire que tout vecteur vérifie la proposition précédente est équivalent à dire que tous les vecteurs de la base vérifient la proposition. Si l'on décompose la solution \bar{u} dans la base des e_i intérieurs, en composantes u_i, on obtient :

 \forall j \in [1,...,n] \sum_{i =1}^n { u_i a( e_i, e_j ) } = \mathcal{L}(e_j)

L'idée est que quand le maillage se resserre et que le nombre de fonctions de base n tend vers l'infini (et que l'espace engendré par cette base  V_n^0 croit vers V0), les solutions un devront converger vers la solution u de l'équation aux dérivées partielles de départ.

Éventuelle deuxième discrétisation

Dans certains problèmes physiques, il peut être intéressant de discrétiser une deuxième fois. Cette seconde discrétisation n'est pas nécessaire pour la méthode des éléments finis. Souvent on a comme expression de \mathcal{L} :  \mathcal{L}(e_j)=\int_\Omega f\cdot e_j

On projette alors f sur la base b. On obtient :  f_n = \sum_{k=1}^n f_k e_k et on approche \mathcal{L}(e_j)= \int_\Omega f\cdot e_j par  \sum_{k=1}^n f_k \int_\Omega e_k e_j,\; j=1,\dots,n.

Le problème est d'obtenir ensuite une projection f_n acceptable sachant qu'il n'y a pas nécessairement de produit scalaire associé à la base qui permette de projeter de façon efficace. Dans les deux exemples de bases donnés plus haut, cette projection est aisée. Dans le cas des éléments finis de Lagrange, la projection sur la fonction ei est donnée par la valeur en xi ; dans le cas des éléments d'Hermite, c'est la valeur de la fonction ainsi que de sa dérivée qui permettent d'obtenir la projection. Pour d'autres bases, la projection peut être plus compliquée.

Problème sous forme matricielle

Si l'on note:

  • la matrice A ayant pour composantes les a( e_i, e_j ),
  • le vecteur U ayant pour composantes les u_i qui sont les coordonnées de la solution approchée sur la base b
  • le vecteur B ayant pour composantes les \mathcal{L}( e_j )

alors ce problème revient à résoudre l'équation linéaire de n équations à n inconnues : \mathbf{A} U = B

La matrice A est appelée matrice de rigidité par analogie avec certains problèmes de mécanique des solides. A est par construction symétrique, et puisque a est coercive, alors A est symétrique, définie positive donc inversible. On obtient donc l'existence et l'unicité de U = A^{-1} B. grâce aux coordonnées de \bar{u} sur la base b on peut alors construire la solution approchée \bar{u}. Quand le maillage se resserre cette solution approchée va tendre vers la vraie solution de l'équation aux dérivées partielles de départ.

Pour le cas avec une deuxième discrétisation de \mathcal{L}( e_j ) on obtient :  \mathbf{A}U=\mathrm{M}f où M est appelée la matrice de masse et contient les \int_\Omega e_i \cdot e_j . f est un vecteur contenant les coordonnées de f dans la base. La méthode est alors la même qu'avec une seule discrétisation puisque A vérifie les mêmes propriétés. Cette méthode peut parfois être préférée quand on peut obtenir de façon simple la projection de f sur la base et la matrice M.

Algorithme

La méthode des éléments finis doit être conduite ainsi

  1. On calcule la matrice de rigidité A
  2. On détermine le membre de droite, en calculant les termes \mathcal{L}(e_j) ou alors par l'intermédiaire de la matrice de masse.
  3. On résout le problème AU=B ou le problème AU=Mf suivant le niveau de discrétisation choisi. U est alors donné par U = A^{-1} B. Selon la base qui a été choisie et selon les données du problème, il faut choisir la méthode d'inversion la plus efficace pour A. C'est l'étape la plus consommatrice en termes de puissance de calcul, et l'efficacité de la méthode en termes de temps de calcul se joue principalement sur cette étape.
  4. On peut écrire \bar{u} grâce au vecteur U qui contient les coordonnées de \bar{u} sur la base b et obtenir une solution approchée au problème.

Condition de Neumann

La condition qui suit est très différente de celle de Dirichlet. On pose comme condition au bord que la dérivée normale existe sur le bord \bar{\Omega}, et que la condition de Neumann \frac{\partial u}{\partial n}(\overline{\Omega}) = 0 est vérifiée.

Si la fonction est supposée différentiable au bord D_u(X) \cdot n = 0, \forall X \in \overline{\Omega}, voire si elle admet un gradient \langle \nabla u(X) | n \rangle = 0, \forall X \in \overline{\Omega}.

Le résultat fonctionne de la même manière car l'élément clef de la démonstration où intervient l'hypothèse de bord est que l'intégrale \int_{\partial \Omega }{\frac{\partial u}{\partial n} v \,\mathrm{d}\omega}=0 parce que cette fois-ci ce n'est pas la fonction test mais la dérivée normale qui est nulle.

Par la suite, la différence réside surtout dans le choix des vecteurs de base pour la discrétisation : il faut conserver les fonctions tests propres aux nœuds du bord.

Exemples issus de problèmes physiques

Historique

La méthode des éléments finis est apparue avec l'analyse des structures, née vers 1850. Les premières études menées sur la résistance des matériaux dans des conditions de petites déformations, ce qui a permis d'obtenir des systèmes simples résolus « manuellement », notamment par Maxwell, Castigliano, Mohr. La formalisation, et le concept mathématique d'élément fini est apparu bien plus tard, vers 1940 et la définition est posée par Newmark, Hrenikoff, Mc Henry et Courant.

L'arrivée du calcul numérique et de méthodes de résolution performantes par ordinateur a permis de populariser la méthode.

Exemple de problème discret : un réseau électrique

Problème discret

Équation locale du composant e : \left. \begin{matrix} I_i^e = {1 \over R^e} (U_i^e - U_j^e) \\ I_j^e = {1 \over R^e} (U_j^e - U_i^e) \end{matrix} \right\} \Longrightarrow \begin{Bmatrix} I_i^e \\ I_j^e \end{Bmatrix} = {1 \over R^e} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}\begin{Bmatrix} U_i^e \\ U_j^e \end{Bmatrix} \, soit \left\{ I^e \right\} = \left[ K^e \right] \left\{ U^e \right\}. \,

On écrit :

  • La continuité des potentiels en chaque connexion  i \,
  • L'équilibre des courants à chaque connexion
  • L'adjonction des courants externes  P^i \,

On obtient l'équation globale du système assemblé : P^i = \sum_{j=1}^m \sum_{e=1}^m K_{ij}^e U_j^e \, soit : \left\{ P \right\} = \left[ K \right] \left\{ U \right\} \,

Définition d'un élément fini

Milieu continu
Milieu discrétisé


En calcul de structures, un élément fini est caractérisé par deux matrices :

  • La matrice de raideur \left[ K \right] \,
  • La matrice de masse \left[ M \right] \,

Application de la méthode des éléments finis à la mécanique

Définitions et notations

On cherche ici à déterminer le vecteur des déplacements \left\{ U \right\}. \, C'est un vecteur dont chaque composante est également appelée degré de liberté (ddl)

  • 3 ddl de translation :  U_x,U_y,U_z \,
  • 3 ddl de rotation :  \theta_x,\theta_y,\theta_z \,

On écrit alors le tenseur des déformations \left[ \epsilon \right] \,, qui modélise la façon dont le matériau va se déformer par rapport à sa position initiale.

Sous l'hypothèse des petites déformations, on a  (A'B')^2 = \left(dx + {\partial u \over \partial x} \ dx\right)^2 + \left({\partial v \over \partial x} \ dx\right)^2 \,

 \epsilon_x = {{A'B'-AB} \over AB} \, Comme  AB=dx \,, on a  A'B'=(1+\epsilon_x)\, dx \,  \Rightarrow 2\epsilon_x+\epsilon_x^2 = 2{\partial u \over \partial x}+\left({\partial u \over \partial x}\right)^2+\left({\partial v \over \partial x}\right)^2 \, Sous l'hypothèse des petites déformations, on néglige les termes d'ordre 2 :  \epsilon_x = {\partial u \over \partial x} \,

Remarque :  \epsilon \, est sans dimension.

On écrit ensuite le tenseur des contraintes) \left[ \sigma \right] \,, qui représente les efforts internes qui s'appliquent dans la structure.

  • \sigma_{11},\sigma_{22},\sigma_{33} \, : contraintes normales
  • \sigma_{12},\sigma_{13},\sigma_{23} \, : contraintes de cisaillement

Une contrainte est homogène à une pression (N/m²).

Ces deux tenseurs sont liés par une loi de comportement. On considère la loi de Hooke : \begin{Bmatrix}\sigma_{11} \\ \sigma_{22} \\ \sigma_{33}
 \\ \sigma_{12} \\ \sigma_{23} \\ \sigma_{31}\end{Bmatrix} = {{E(1-\nu)} \over {(1+\nu)(1-2\nu)}} \begin{bmatrix}1 & {\nu \over {1-\nu}} & {\nu \over {1-\nu}} & 0 & 0 & 0 \\ & 1 & {\nu \over {1-\nu}} & 0 & 0 & 0 \\ & & 1 & 0 & 0 & 0 \\ & & & {(1-2\nu) \over {2(1-\nu)}} & 0 & 0 \\ & \mbox{(sym)} & & & {(1-2\nu) \over {2(1-\nu)}} & 0 \\ & & & & & {(1-2\nu) \over {2(1-\nu)}}\end{bmatrix} \begin{Bmatrix}\epsilon_{11} \\ \epsilon_{22} \\ \epsilon_{33} \\ \gamma_{12} \\ \gamma_{23} \\ \gamma_{31}\end{Bmatrix}\, avec

  • E \,: module de Young (N/m²)
  • \nu \, : coefficient de Poisson (sans dimension)

On utilise parfois le module de cisaillement : G = {E \over {2(1+\nu)}} \,

Pour un matériau isotrope, il n'y a que deux paramètres indépendants. Il y en a 6 pour un matériau isotrope transverse, 9 pour un matériau orthotrope et 21 pour un matériau anisotrope.

En notation matricielle, on écrit : \left\{ \sigma \right\} = \left[D\right]\left\{ \epsilon \right\} \, \left[D\right] \, est appelée matrice d'élasticité du matériau.

L'énergie de déformation W s'écrit W = {1 \over 2}\int_{v}\sigma.\epsilon.\,\mathrm{d}v \,

Travail d'une force
C'est le produit de la force par le déplacement de son point d'application : \tau_F = {1 \over 2}F.U_F \,
Moment
C'est la mesure d'une action entraînant une rotation, proportionnelle à la force et à la longueur du bras de levier entre le centre de rotation du solide et le point d'application de la force.

Équations fondamentales

Équations d'équilibre local
{\partial \sigma_x \over \partial x} + {\partial \sigma_{xy} \over \partial y} + {\partial \sigma_{xz} \over \partial z} + F_x = 0\,
{\partial \sigma_y \over \partial y} + {\partial \sigma_{xy} \over \partial x} + {\partial \sigma_{yz} \over \partial z} + F_y = 0\,
{\partial \sigma_z \over \partial z} + {\partial \sigma_{zx} \over \partial x} + {\partial \sigma_{zy} \over \partial y} + F_z = 0\,
Relations déformations-déplacements
{\epsilon_x= {\partial u \over  \partial x}} \qquad {\epsilon_y= {\partial v \over  \partial y}} \qquad {\epsilon_z= {\partial w \over  \partial z}} \,
\gamma_{xz}={\partial w \over \partial x} + {\partial u \over \partial z} \,
\gamma_{xy}={\partial u \over \partial y} + {\partial v \over \partial x} \,
\gamma_{yz}={\partial w \over \partial y} + {\partial v \over \partial z} \,
Symboliquement, on écrit \{\epsilon\} = [S]\{U\} \,
Signification du coefficient de Poisson

Si on applique au barreau une contrainte \sigma_x \,, on observe un rétrécissement dans la direction y correspondant à une déformation -\nu \sigma_x \over E \,

Exemple de formulation : barre en traction

On suppose que le déplacement en tout point de la barre est donné par un polynôme du 1er degré :  u(x) = a_1 + a_2 x \,

On a  u(0) = u_1 \, et  u(L) = u_2 \,
d'où  u(x) = u_1\left(1-{x \over L} \right) + {x \over L} u_2 \,

qu'on écrit symboliquement :  u(x) = [N(x)]\{U\} \, avec \left|\begin{matrix}[N(x)] & = & \left[ \left(1-{x \over L}\right) ; \left({x \over L}\right) \right] \\ 
\{U\} & = & \begin{Bmatrix}u_1 \\ u_2 \end{Bmatrix} \end{matrix}\right. \,
On en déduit :

 \epsilon = \left[-{1 \over L} ; {1 \over L} \right]\begin{Bmatrix}u_1 \\ u_2\end{Bmatrix} = [B]\{U\} \,
 \sigma = E\left[-{1 \over L} ; {1 \over L} \right]\begin{Bmatrix}u_1 \\ u_2\end{Bmatrix} = [D]\epsilon \,

D'autre part, on a par définition :

 \begin{Bmatrix}F_1 \\ F_2 \end{Bmatrix} = \begin{Bmatrix}-\sigma S \\ \sigma S \end{Bmatrix} \, où S est l'aire de la section de la barre.

On pose :

 \{F\} = \begin{Bmatrix}F_1 \\ F_2 \end{Bmatrix} = [A]\sigma \qquad \mbox {avec} \qquad \{A\} = S \begin{Bmatrix}-1 \\ 1 \end{Bmatrix} \,

On obtient finalement :

 \{F\} = [A][D][B]\{U\} \,

Soit une relation du type :

 \{F\} = [K]\{U\} \qquad \mbox {avec} \qquad [K] = [A][D][B] \,

En explicitant :

 [K] = {ES \over L} \begin{bmatrix}1 & -1 \\ -1 & 1 \end{bmatrix} \,

On voit que la matrice de rigidité se calcule comme le produit de trois matrices :

 [B] \, : Transformation des déplacements aux déformations
 [D] \, : Matrice d'élasticité du matériau
 [A] \, : Transformation des contraintes en forces

Formulation générale (méthode directe)

La démarche est la suivante :

  • On exprime le déplacement \{u(x)\}\, en tout point de l'élément en fonction des déplacements aux nœuds \{u(x)\} = [N(x)]\{U\}\,
  • On exprime les déformations en fonction des déplacements \{\epsilon(x)\}=[S]\{u(x)\},\, d'où
\{\epsilon(x)\}=[S]\{N(x)\}\{U\}=[B(x)]\{U\}\,

On écrit la loi de comportement du matériau qui relie les contraintes aux déformations :

\{\sigma(x)\}=[D]\{\epsilon(x)\}\,

On écrit que le travail des forces externes appliquées à la structure pour un déplacement virtuel \delta U \, est égal au travail interne des contraintes pour ce même déplacement :

\{\delta U\}^T\{F\}=\int_V \{\delta \epsilon\}^T\{\sigma\} \,\mathrm{d}v \,

En explicitant, on a :

\{\delta U\}^T\{F\}=\{\delta U\}^T \left (\int_V [B]^T[D][B] \,\mathrm{d}v \right ) \{U\} \,

Comme cette relation est vraie pour tout déplacement virtuel, on en déduit :

\{F\}=[K]\{U\} \,

avec [K] \, sous sa forme plus générale : [K]=\int_V [B]^T[D][B] \,\mathrm{d}v \,

Remarques
  • La relation ci-dessus montre que [K] \, est symétrique.
  • Le terme courant K_{ij} \, de la matrice correspond à la force qui s'exerce sur le nœud j \, lorsqu'on impose un déplacement unitaire du nœud i \,.

La symétrie de [K] \, qui s'écrit K_{ij} = K_{ji} \, correspond mécaniquement au théorème de réciprocité de Maxwell-Betti.

Qu'est-ce qui va différencier les différents types d'éléments finis ?

  • Le choix des fonctions N(x)
  • La nature de l'opérateur S (reliant déformations et déplacements) qui dépend du type de théorie élastique utilisée :
    • théorie des poutres
    • théorie des plaques
    • théorie de la contrainte ou de la déformation plane
    • théorie des coques
    • théorie des corps de révolution
    • théorie de l'élasticité 3D
Remarques

Le processus de formulation d'un élément fini décrit ici est celui de la méthode directe (dite aussi méthode des déplacements). Il existe d'autres approches :

  • La méthode des résidus pondérés
  • L'application du principe des travaux virtuels ou des puissances virtuelles
  • La minimisation de l'énergie potentielle

Toutes ces approches sont équivalentes et aboutissent à la construction de la même matrice de rigidité.

Éléments finis en contrainte

Au lieu de rechercher une solution approchée en déplacement, on peut aussi rechercher la solution approchée en contrainte.

Dans le cas de la mécanique, l'application du principe des puissances virtuelles donne de manière non triviale les théorèmes énergétiques. On peut aboutir au même résultat en quelques lignes en écrivant l'erreur en relation de comportement.

L'approche en contrainte consiste à rechercher dans l'espace des champs de contraintes admissibles celui qui réalise le minimum de l'énergie complémentaire.

Cette approche est plus précise que l'approche en déplacement mais elle est peu développée du fait de la difficulté que l'on a à générer des champs de contraintes de divergence donnée.

Étude des fonctions N(x)

  • Dans le cas général de l'élasticité tridimensionnelle, ce sont en fait des fonctions de x, y, z.
  • Les fonctions les plus couramment utilisées sont des polynômes.
    • Polynôme de degré 1 : élément linéaire (2 nœuds par arête)
    • Polynôme de degré 2 : élément parabolique (3 nœuds par arête)
    • Polynôme de degré 3 : élément cubique (4 nœuds par arête)

Les fonctions N(x) sont appelées fonctions de forme ou fonctions d'interpolation de l'élément.

Les éléments isoparamétriques

[K]=\int_V [B]^T[D][B] \,\mathrm{d}v \,

Il arrive alors un dilemme : soit on construit [K] \, pour un certain nombre d'éléments de forme et de géométrie figée, et il est alors nécessaire, pour mailler une structure complexe, d'utiliser un grand nombre d'éléments, soit on utilise des éléments à géométries variables, et il faut reconstruire [K] \, à chaque fois.

Une solution courante est alors d'utiliser des fonctions d'interpolation pour décrire non seulement le champ de déplacement de l'élément mais également sa géométrie, tout en travaillant en coordonnées locales.

Interpolation de la géométrie

\begin{matrix}x &=&\ &\tilde N_1(x)x_1 + \tilde N_2(x)x_2 \\
\ &\ &+&\tilde N_3(x)x_3 + \tilde N_4(x)x_4\end{matrix} \,

Idem pour les autres coordonnées.

Coordonnées locales (cas 2D)

Élément isoparamétrique

Un élément est dit isoparamétrique si on prend les mêmes fonctions d'interpolation pour le déplacement et la géométrie.
\tilde N(x) = N(x) \,

Autres classes d'éléments

Évaluation de K

La forme générale s'écrit :
[K]=\int_x \int_y \int_z G(x,y,z) \,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z \,
On passe en variables locales \xi,\eta,\zeta \,
On a \,\mathrm{d}x\,\mathrm{d}y\,\mathrm{d}z = \det (J)\,\mathrm{d}\xi \,\mathrm{d}\eta \,\mathrm{d}\zeta \,

J \, s'appelle la matrice jacobienne.

On est alors amené à calculer des intégrales du type :
\int_{-1}^{+1} \int_{-1}^{+1} \int_{-1}^{+1} G(\xi ,\eta ,\zeta) \det J\ \,\mathrm{d}\xi \,\mathrm{d}\eta \,\mathrm{d}\zeta \,

Bénéfice de l'approche

On s'est ramené à un domaine d'intégration simple et invariant pour lequel on peut appliquer les formules de quadrature de Gauss :
\int_{-1}^{+1} f(x) \,\mathrm{d}x = \sum_{k=1}^n H_k f(a_k) \,      les a_k \, et H_k \,étant tabulés.
Les a_k \, sont appelés points d'intégration de l'élément ou encore points de Gauss de l'élément.

Cas particulier : les éléments axisymétriques

Décomposition en série de Fourier :

u(r,\theta,z) = \sum_n u_n^s(r,z)\cos (n\theta) + \sum_n u_n^s(r,z)\sin (n\theta) \,

L'axisymétrie correspond à la restriction n=0 \, de cette décomposition.

u(r,\theta,z) = u_n^s (r,z) \,
Remarque

Pour utiliser ce type d'élément, le problème doit être globalement axisymétrique :

  • la géométrie
  • les conditions limites
  • le chargement

Processus de calcul (cas statique)

  1. Maillage
  2. Construction de la matrice de raideur de chaque élément [K^e] \,
  3. Assemblage de la matrice globale [K] \,
  4. Construction du vecteur chargement \{F\} \,
  5. Élimination de certains ddl (si besoin)
  6. Résolution : \{U\} = [K^{-1}]\{F\} \,
  7. Calcul des quantités dérivées de \{U\} \,
    \{\epsilon^e\} = [B^e] \{U^e\} \,
    \{\sigma^e\} = [D^e] \{\epsilon^e\} \,
    W^e = {1 \over 2} \{\epsilon^e\}^T \{\sigma^e\} \,

Logiciels d'éléments finis

Quelques exemples de logiciels utilisant la méthode des éléments finis en mécanique des structures :

  • ABAQUS : logiciel pluridisciplinaire développé par la société Dassault Systèmes
  • ADVANCE DESIGN : logiciel français développé par GRAITEC ( remplace EFFEL)
  • ANSYS : logiciel pluridisciplinaire développé par ANSYS
  • CAST3M : logiciel pluridisciplinaire français développé par le CEA (gratuit pour l'enseignement et la recherche)
  • ASTER : logiciel pluridisciplinaire libre français développé par EDF
  • COMSOL Multiphysics : logiciel pluridisciplinaire développé par Comsol,
  • CosmosWorks : Ancienne version de SolidWorks (Dassault Systèmes)
  • Dytran : logiciel américain développé par MSC.Software
  • EFFEL : logiciel français développé par GRAITEC
  • EuroPlexus : logiciel français
  • Flux2D/3D : logiciel 2D&3D éléments finis français (développé en collaboration avec le GE2Lab) permettant le calcul des états magnétiques, électriques ou thermiques en régimes permanents, harmoniques et transitoires, avec des fonctionnalités d'analyse multi-paramétrique étendue, les couplages circuit et cinématique.
  • FreeFem++ : logiciel pluridisciplinaire libre.
  • Feel++ : logiciel pluridisciplinaire libre.
  • Getfem++ : logiciel pluridisciplinaire libre (avec dominante mécanique des structures).
  • ICAB : logiciel de calcul pour poutres et coques avec vérifications réglementaires selon Document technique unifié DTU France (NV65, CM66, AL76, CB71...), Eurocodes, AISC américain
  • IMPETUS Afea Solver : logiciel de calcul de grandes déformations utilisant des éléments finis isogéométriques d'ordre 3.
  • JMAG : logiciel japonais (distribué en Europe par Powersys) permet un couplage entre les analyses électromagnétiques et structurelles.
  • LS-DYNA : logiciel de dynamique rapide américain
  • Marc : logiciel de grandes déformations américain développé par MSC.Software
  • Morfeo : logiciel belge
  • MSC.Nastran : logiciel Pluridisciplinaire américain développé par MSC.Software
  • NX.Nastran : logiciel développé par Siemens
  • PAM-CRASH : logiciel de dynamique rapide (crash) développé par la société française ESI
  • PAM-STAMP : logiciel de simulation d'emboutissage et de mise en forme de tôle, développé par la société française ESI
  • PERMAS : logiciel pluridisciplinaire développé par la société allemande INTES GmbH, basée à Stuttgart.
  • Radioss : logiciel pluridisciplinaire développé par la société américaine Altair
  • ROBOT MILLENIUM : logiciel français développé par ROBOBAT pour le calcul de structures de type génie civil et bâtiment, il est intégré ensuite par AutoDesk
  • RFEM : logiciel de calcul de structure, développé par la société allemande Dlubal Software.
  • SAMCEF : logiciel pluridisciplinaire belge, développé par SAMTECH, propriété de Siemens Automotive
  • Structurix : Programme gratuit de calcul de résistance des matériaux. Il permet de résoudre des problèmes 2D et 3D avec des éléments : barres et poutres. Des problèmes 2D avec des éléments triangulaires. Des problèmes de flexion des plaques à l'aide d'éléments rectangulaires.
  • Z-set : logiciel pluridisciplinaire développé par MINES ParisTech, l'ONERA et NorthWest Numerics, spécialisé dans les comportements non linéaires des matériaux
  • SYSWELD : logiciel de Thermo-mécano-métallurgie français basé sur SYSTUS développé par la société ESI

Notes et références

  1. Pour certains cas limites lorsqu'il n'y a pas suffisamment de régularité et que l'intégration par parties n'a pas de sens on utilise aussi une formulation adaptée du théorème de Stokes mais il faut alors définir la dérivée faible puis la divergence faible, ce qui entraine des complications supplémentaires.
  2. (en) Advanced Finite Element Methods, University of Colorado at Boulder (lire en ligne [PDF]), « The Linear Tetrahedron »
  3. (en) Advanced Finite Element Methods, University of Colorado at Boulder (lire en ligne [PDF]), « The Quadratic Tetrahedron »
  4. (en) Advanced Finite Element Methods, University of Colorado at Boulder (lire en ligne [PDF]), « Hexahedron Elements »

Voir aussi

Articles connexes

  • Calcul de structures et modélisation
  • NAFEMS
  • Méthode des différences finies
  • Méthode des volumes finis
  • Méthode des éléments finis de frontière
  • Théorème de Lax-Milgram
  • Méthode des éléments finis étendus
  • Analyse isogéométrique

Bibliographie

  • G. Allaire and A. Craig: Numerical Analysis and Optimization:An Introduction to Mathematical Modelling and Numerical Simulation
  • K. J. Bathe : Numerical methods in finite element analysis, Prentice-Hall (1976) (ISBN 0136271901)
  • J. Chaskalovic, Méthode des éléments finis pour les sciences de l'ingénieur (2004), Éd. Lavoisier. (ISBN 2-7430-0708-7)
  • Chu-Kia Wang : Matrix methods of structural analysis, International Textbook Co ; 2e édition (1970) (ISBN 0700222677)
  • J. C. Cuillière, Introduction à la méthode des éléments finis (2011), Éd. Dunod. (ISBN 2100564382)
  • R. H. Gallagher : Introduction aux éléments finis, Pluralis (1977)
  • N. Willems : Matrix analysis for structural engineers, Prentice-Hall (1968) (ISBN 0135654998)
  • O. C. Zienkiewicz, R. L. Taylor, J. Z. Zhu : The Finite Element Method: Its Basis and Fundamentals, Butterworth-Heinemann ; 6e édition (21 mars 2005) (ISBN 0750663200)

Liens externes

  • (fr) RDM Le Mans v.6, logiciel téléchargeable
  • (en) La méthode des éléments finis en génie biomedicale et en biomécanique
  • (fr) Méthode des éléments finis, Hervé Oudin, École Centrale de Nantes
  • (fr) La Méthode des Éléments Finis: Vulgarisation des aspects mathématiques, Illustration des capacités de la méthode, Vincent Manet
  • (fr) Méthodes numériques appliquées à la conception par éléments finis, David Dureisseix, Université Montpellier II - Sciences et Techniques du Languedoc
  • (fr) Formulations mathématiques et résolution numérique en mécanique, Christian Wielgosz, Bernard Peseux, Yves Lecointe, Université de Nantes
  • (fr) Design Solutions for Electrical Engineering,Flux pionnier dans le développement des formulations EF
  • Portail de l’analyse
  • Portail de la physique
  • Portail du génie mécanique
This article is issued from Wikipédia - version of the Friday, September 25, 2015. The text is available under the Creative Commons Attribution/Share Alike but additional terms may apply for the media files.
Contents Listing Alphabetical by Author:
A B C D E F G H I J K L M N O P Q R S T U V W X Y Z Unknown Other

Contents Listing Alphabetical by Title:
# A B C D E F G H I J K L M N O P Q R S T U V W Y Z Other

Medical Encyclopedia

Browse by first letter of topic:


A-Ag Ah-Ap Aq-Az B-Bk Bl-Bz C-Cg Ch-Co
Cp-Cz D-Di Dj-Dz E-Ep Eq-Ez F G
H-Hf Hg-Hz I-In Io-Iz J K L-Ln
Lo-Lz M-Mf Mg-Mz N O P-Pl Pm-Pz
Q R S-Sh Si-Sp Sq-Sz T-Tn To-Tz
U V W X Y Z 0-9

Biblioteca - SPANISH

Biblioteca Solidaria - SPANISH

Bugzilla

Ebooks Gratuits

Encyclopaedia Britannica 1911 - PDF

Project Gutenberg: DVD-ROM 2007

Project Gutenberg ENGLISH Selection

Project Gutenberg SPANISH Selection

Standard E-books

Wikipedia Articles Indexes

Wikipedia for Schools - ENGLISH

Wikipedia for Schools - FRENCH

Wikipedia for Schools - SPANISH

Wikipedia for Schools - PORTUGUESE

Wikipedia 2016 - FRENCH

Wikipedia HTML - CATALAN

Wikipedia Picture of the Year 2006

Wikipedia Picture of the Year 2007

Wikipedia Picture of the Year 2008

Wikipedia Picture of the Year 2009

Wikipedia Picture of the Year 2010

Wikipedia Picture of the Year 2011