Eléments finis
(brève introduction)

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 ₀¹ [u_,xx(x)v(x)]dx + Int ₀¹ [c(x)u(x)v(x)]dx = Int ₀¹ [f(x)v(x)]dx

Une intégration par parties du premier terme donne :

Int ₀¹ [u_,x(x)v_,x(x)]dx - u_,x(1)v(1) + u_,x(0)v(0) + Int ₀¹ [c(x)u(x)v(x)]dx = Int ₀¹ [f(x)v(x)]dx

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

Int ₀¹ [u_,x(x)v_,x(x)]dx + Int ₀¹ [c(x)u(x)v(x)]dx = Int ₀¹ [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 w₁, .. ,w_N N fonctions linéairement indépendantes de V. On note V_h le sous-espace vectoriel de V engendré par les w_i. On peut approximer le problème variationnel par le suivant : trouver uh appartenant à Vh telle que

Int ₀¹ [u_h,x(x)v_h,x(x)]dx + Int ₀¹ [c(x)u_h(x)v_h(x)]dx = Int ₀¹ [f(x)v_h(x)]dx

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

u_h(x) = somme_i u_iw_i(x)

où les u_i sont N réels à déterminer. En prenant v_h = w_j, on cherche alors les u_i tels que :

somme_i u_i (Int ₀¹ [w_i,xw_j,x]dx + Int ₀¹ [c(x)w_i(x)w_j(x)].dx) = Int ₀¹ [f(x)w_j(x)]dx

Si on note A la matrice NxN de coefficients

A_ij = (Int ₀¹ [w_i,xw_j,x]dx + Int ₀¹ [c(x)w_i(x)w_j(x)].dx)

et f le vecteur de composantes Int ₀¹ [f(x)w_j(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 w_i de manière à ce que la matrice A soit creuse et que la solution u_h 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 x_i = ih, où h = 1/(N+1). On définit les fonctions w_i(x) :

(x-x_i-1)/(x_i-x_i-1) si x_i-1 < x < x_i
w_i(x) = (x-x_i+1)/(x_i-x_i+1) si x_i < x < x_i+1
0 sinon

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

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

2/h si i=j
Int ₀¹ [w_i,xw_j,x]dx = -1/h si |i-j| =1
0 sinon

Le calcul de Int ₀¹ [c(x)w_i(x)w_j(x)] et de Int ₀¹ [f(x)w_j(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 :

-Lap (u) = q sur un domaine V
u = 0 sur une partie S₁ de sa frontière S
u_,n = g sur le complément S₂

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 S₁). En effet, si l'on part d'un problème de Laplace Lap(u) = 0, u=f sur S₁ et u_,n = g sur S₂, et que l'on dispose d'une solution u₀ de ce problème telle que u₀ = f sur S₁, on peut, en posant u₁ = u - u₀, retomber sur le problème de Poisson ci-dessus avec q = Lap(u₀) et en remplaçant g par g - u_0,n. Passer de f sur S₁ à u₀ 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 :

Int_V [- v Lap (u)]dV = Int_V [(grad v).(grad u)]dV - Int_S [u_,n v]dS

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

a(u,v) = Int_V [q v]dV + Int_S2 [g v]dS

où

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

a est une forme bilinéaire symétrique. Et L(v) = Int_V [q v]dV + Int_S2 [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 A_i (i=1,..., N_s) 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 u_p, u_q et u_r aux trois sommets A_p(x_p,y_p), A_q(x_q,y_q) et A_r(x_r,y_r). A, B et C sont les solutions du système linéaire 3x3 : Ax_p + By_p + C = u_p

Ax_q + By_q + C = u_q
Ax_r + By_r + C = u_r

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 l₁, l₂ et l₃ 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 w_i telle que w_i = 1 en A_i et w_i = 0 en A_j si j différent de i. On a alors

u(x,y) = somme_i u_iw_i(x,y), où u_i = u(A_i), i = 1, ..., N_s

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

somme_j a(w_j,w_i)u_j = L(w_i), i = 1, ...., N

où N < N_s, les inconnues u_i pour i compris entre N et Ns étant imposées par les conditions aux limites ... La matrice A formée par les coefficients a(w_i,w_j) 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 w_i est constitué des seuls triangles admettant A_i comme sommet.