Elements finis


Revenons sur un problème aux limites unidimensionnel pour introduire de manière simplifiée la méhode des éléments finis. Le moment fléchissant u(x) d'une poutre de longueur unité, en appui en x = 0 et x = 1, soumise à une force axiale P et à une densité de charge transverse f(x), est solution de l'équation

-u,xx + c(x) u(x) = f(x)
u(0) = u(1) = 0

où c(x) = P/EI(x), avec E le module d'Young et I le moment d'inertie au point x.

On peut comme précédemment résoudre ce problème en choisissant un schéma aux différences finies. La méthode de Galerkin choisit une autre voie : multiplions le problème par une fonction v continûment dérivable sur [0,1]. En intégrant le résultat sur ce même intervalle, on obtient :

- Int 01 [u,xx(x)v(x)]dx + Int 01 [c(x)u(x)v(x)]dx = Int 01 [f(x)v(x)]dx

Une intégration par parties du premier terme donne :

Int 01 [u,x(x)v,x(x)]dx - u,x(1)v(1) + u,x(0)v(0) + Int 01 [c(x)u(x)v(x)]dx = Int 01 [f(x)v(x)]dx

Si l'on impose à la fonction v d'être nulle en 0 et 1, on en déduit :

Int 01 [u,x(x)v,x(x)]dx + Int 01 [c(x)u(x)v(x)]dx = Int 01 [f(x)v(x)]dx

Notons V l'ensemble des fonctions g continues, dont la dérivée première g,x est continue par morceaux (i.e continue sauf éventuellement en un nombre fini de points où g,x possèderait des limites à droite et à gauche). On peut montrer que V est un espace vectoriel. L'objectif est de chercher les fonctions u appartenant à V qui vérifient la formulation variationnelle que l'on vient d'obtenir pour tout v appartenant à V.

Les solutions du problème variationnel sont à priori moins régulières, puisque le problème différentiel initial contenait des dérivées secondes, alors que le problème variationnel ne fait plus intervenir que des dérivées premières. On peut cependant montrer que la solution de l'un des problèmes est solution de l'autre, et réciproquement, si c(x) 0.

Soient w1, .. ,wN N fonctions linéairement indépendantes de V. On note Vh le sous-espace vectoriel de V engendré par les wi. On peut approximer le problème variationnel par le suivant : trouver uh appartenant à Vh telle que

Int 01 [uh,x(x)vh,x(x)]dx + Int 01 [c(x)uh(x)vh(x)]dx = Int 01 [f(x)vh(x)]dx

Ceci constitue l'approximation de Galerkin du problème faible. On peut écrire

uh(x) = sommei uiwi(x)

où les ui sont N réels à déterminer. En prenant vh = wj, on cherche alors les ui tels que :

sommei ui (Int 01 [wi,xwj,x]dx + Int 01 [c(x)wi(x)wj(x)].dx) = Int 01 [f(x)wj(x)]dx

Si on note A la matrice NxN de coefficients

Aij = (Int 01 [wi,xwj,x]dx + Int 01 [c(x)wi(x)wj(x)].dx)

et f le vecteur de composantes Int 01 [f(x)wj(x)]dx, cela revient à chercher u tel que

A.u = f

Comme dans la méthode des différences finies, il convient ensuite de résoudre le système linéaire obtenu.

La méthode des éléments finis consiste à choisir les fonctions wi de manière à ce que la matrice A soit creuse et que la solution uh du problème approché converge vers la solution du problème faible.

Pour le problème ci-dessus, divisons [0,1] en N+1 parties de même longueur, et posons xi = ih, où h = 1/(N+1). On définit les fonctions wi(x) :

wi est donc une fonction affine par morceaux, nulle partout sauf sur l'intervalle [xi-1,xi+1] où elle vaut 0 en xi-1, 1 en xi et 0 en xi+1.

La résolution du système linéaire passe par le calcul de la matrice A. On a :

Le calcul de Int 01 [c(x)wi(x)wj(x)] et de Int 01 [f(x)wj(x)]dx peut se faire par intégration numérique, par exemple par la méthode des trapèzes. Le système obtenu dans ce cas est identique à celui obtenu par la méthode des différences finies. La méthode des éléments finis offre cependant davantage de souplesse au niveau :
- de la modélisation des conditions aux limites
- de la répartition non-uniforme des noeuds de discrétisation afin de mieux " suivre " les variations importantes de la solution
- de la possibilité d'introduire des fonctions de degré supérieur à 1, afin d'augmenter la précision de la méthode.
 

Problème elliptique

Soit le problème de Poisson déjà introduit, qui modélise par exemple le déplacement transverse u d'une membrane tendue et soumise à une densité de force verticale proportionnelle à q :

Le second membre q de l'équation de Poisson vient de l'intérêt de disposer d'une condition de Dirichlet homogène (i.e d'une donnée nulle sur S1). En effet, si l'on part d'un problème de Laplace Lap(u) = 0, u=f sur S1 et u,n = g sur S2, et que l'on dispose d'une solution u0 de ce problème telle que u0 = f sur S1, on peut, en posant u1 = u - u0, retomber sur le problème de Poisson ci-dessus avec q = Lap(u0) et en remplaçant g par g - u0,n. Passer de f sur S1 à u0 dans V s'appelle " faire un relèvement ". On verra ultérieurement quel est l'intérêt de disposer d'une condition de Dirichlet homogène.

Multiplions le problème ci-dessus par une fonction v et intégrons le résultat obtenu sur le domaine V. Avec la première identité de Green, qui dit que :

IntV [- v Lap (u)]dV = IntV [(grad v).(grad u)]dV - IntS [u,n v]dS

on obtient le problème variationnel : trouver u tel que, pour tout v,

a(u,v) = IntV [q v]dV + IntS2 [g v]dS

a(u,v) = IntV [grad u. grad v]dV

a est une forme bilinéaire symétrique. Et L(v) = IntV [q v]dV + IntS2 [g v]dS est également une forme linéaire. On peut montrer que la solution de ce problème minimise la fonctionnelle J(v) = 1/2 a(v,v) - L(v). L'introduction de cette fonctionnelle n'est pas à proprement parler nécessaire, mais elle correspond parfois directement à la formulation du problème physique, comme par exemple en élastostatique linéaire.

Supposons que notre domaine V soit bidimensionnel. On peut alors définir une triangulation de V, et noter Ai (i=1,..., Ns) les sommets des triangles obtenus. On va chercher une solution polynomiale par morceaux, c'est-à-dire définie par u(x) = Ax + By + C sur chaque triangle, et déterminer A, B et C en fonction des valeurs up, uq et ur aux trois sommets Ap(xp,yp), Aq(xq,yq) et Ar(xr,yr). A, B et C sont les solutions du système linéaire 3x3 : Axp + Byp + C = up

Axq + Byq + C = uq
Axr + Byr + C = ur

Si aucun triangle n'est aplati, le déterminant de ce système est non-nul et l'on peut bien définir une fonction du premier degré par morceaux de façon unique par ses valeurs aux trois sommets d'un triangle. On peut définir trois fonctions l1, l2 et l3 de x et y affines et telles que chacune vaille 1 sur l'un des sommets et 0 sur les deux autres. On peut alors définir ce que l'on nomme les coordonnées barycentriques d'un point M du triangle considéré. En se plaçant sur le domaine entier, on peut de la même manière définir une base canonique wi telle que wi = 1 en Ai et wi = 0 en Aj si j différent de i. On a alors

u(x,y) = sommei uiwi(x,y), où ui = u(Ai), i = 1, ..., Ns

On peut alors remplacer le problème variationnel par le problème approché a(u,wi) = L(wi) qui mène à :

sommej a(wj,wi)uj = L(wi), i = 1, ...., N

où N < Ns, les inconnues ui pour i compris entre N et Ns étant imposées par les conditions aux limites ... La matrice A formée par les coefficients a(wi,wj) est symétrique et définie positive du fait des propriétés de a, et l'on peut ramener l'écriture précédente à celle d'un système Ax = b. A est également très creuse, car le support de chaque fonction wi est constitué des seuls triangles admettant Ai comme sommet.