4.2 Formulation éléments finis

Pour construire une formulation variationnelle discrète à partir de la formulation (5.6), il faut construire une approximation $v^{h}$ de la solution exacte, vérifiant les conditions cinématiques (conditions de Dirichlet):

$$v^{h}(0)=0, \frac{dv^{h}}{dx}(0)=0$$ (4.7)

et telle que l'on puisse calculer le Lagrangien discret associé à (5.5):

$$\mathcal{L}(v^{h})=-\frac{1}{2}\int_{0}^{L}EI_{zz}\left(\fra... ...right)^{2} dx+\left(\int_{0}^{L}qv^{h} dx+Fv^{h}(L)\right)$$ (4.8)

Pour construire cette approximation par éléments finis , on impose les contraintes suivantes:

1. la solution $v^{h}$ doit être continue ainsi que sa dérivée $\frac{dv^{h}}{dx}$ sur tout le domaine $\Omega=[0,L]$
2. sur chaque élément, la solution est un polynôme

La première condition permet le calcul des intégrales dans (5.8), puisque $v^{h}$ est $\mathcal{C}^{1}(\Omega)$ et permet d'imposer les conditions aux limites fortes (5.7). La seconde condition permet de construire une approximation simple. Sur un élément la solution est un polynôme, qui doit vérifier quatre conditions: la continuité de la fonction et de sa dérivée en chacune des deux extrémités de l'élément. C'est donc un polynôme d'Hermite de degré 3. Sur la figure 5.1, on a tracer un exemple d'approximation pour un maillage de 2 éléments.

Figure 5.1: approximation $v^{h}$ sur un maillage de 2 éléments

D'un point de vue mécanique, on cherche une solution approchée $v^{h}$ qui fournit une approximation continue du déplacement $v(x)$ et de l'angle de rotation $\theta(x)$ (dérivée première de $v(x)$), et une approximation par élément (éventuellement discontinue aux noeuds) du moment fléchissant $M_{z}$ (dérivée seconde de $v(x)$) et de la résultante de la contrainte de cisaillement $V_{y}$ (dérivée troisième de $v(x)$).

On construit donc un maillage de $\Omega$ correspondant à un découpage en $ne$ éléments:

$$\Omega=\bigcup_{i=1}^{ne}[x_{i-1},x_{i}]$$

soient $ne+1$ noeuds ${x_{i}}_{i=0..ne}$.

L'approximation utilise deux degré de liberté par noeuds: pour chaque noeud $i$ de coordonnée $x_{i}$ ce sont la valeur nodale de la fonction $v_{i}=v^{h}(x_{i})$ et la valeur nodale de la dérivée $(\frac{dv}{dx})_{i}=\frac{dv^{h}}{dx}(x_{i})$. Cela permet de définir de façon unique l'approximation $v^{h}$, qui possède donc $2ne=2nn$ degrés de liberté (les valeurs au noeud $i=0$ sont imposées par la condition aux limites 5.7). Cette approximation s'écrit comme une combinaison linéaire des valeurs nodales de la fonction et de sa dérivée.

En notant $\Phi_{i}(x)$ les fonctions de base associées aux valeurs nodales de la fonction $v_{i}$ et $\Psi_{i}(x)$ les fonctions de base associées aux valeurs nodales de la dérivée $(\frac{dv}{dx})_{i}$, on écrit:

$$v^{h}(x)=\sum_{i=1}^{nn}v_{i}\Phi_{i}(x)+\sum_{i=1}^{nn}(\frac{dv}{dx})_{i}\Psi_{i}(x)$$

Sur un élément $e_{k}=[x_{k-1},x_{k}]$, cette approximation s'écrit:

$$v^{h}(x)=v_{k-1}\Phi_{k-1}(x)+(\frac{dv}{dx})_{k-1}\Psi_{k-1}(x)+v_{k}\Phi_{k}(x)+(\frac{dv}{dx})_{k}\Psi_{k}(x)$$

C'est l'unique polynôme de degré $\le3$, $p(x)$ qui vérifie:

$$p(x_{k-1})=v_{k-1}, p(x_{k})=v_{k}, \frac{dp}{dx}(x_{... ...rac{dv}{dx})_{k-1}, \frac{dp}{dx}(x_{k})=(\frac{dv}{dx})_{k}$$

C'est un polynôme d'Hermite de degré 3, d'où le nom d'éléments finis hermitiens que l'on donne à cette approximation.

4.2.1 fonctions de forme

Pour calculer l'approximation sur un élément $e_{k}$, on se place sur l'élément de référence $[-1,+1]$ (comme dans le chapitre précédent) en effectuant un changement de variable:

$$\xi=\frac{2x-x_{k}-x_{k-1}}{x_{k}-x_{k-1}}$$ (4.9)

Sur cet élément de référence, on détermine les 4 fonctions de forme associées, qui sont les polynômes de Hermite de degré 3:

1. $N_{1}(\xi)$ associée à la valeur nodale au premier noeud ($\xi=-1$): i.e. telle que
   
   
   $$N_{1}(-1)=1, N_{1}(1)=0, \frac{dN_{1}}{d\xi}(-1)=0, \frac{dN_{1}}{d\xi}(1)=0$$
   
   

2. $N_{2}(\xi)$ associée à la valeur nodale au second noeud ($\xi=+1$): i.e. telle que
   
   
   $$N_{2}(-1)=0, N_{2}(1)=1, \frac{dN_{2}}{d\xi}(-1)=0, \frac{dN_{2}}{d\xi}(1)=0$$
   
   

3. $H_{1}(\xi)$ associée à la valeur de la dérivée au premier noeud ($\xi=-1$): i.e. telle que
   
   
   $$H_{1}(-1)=0, H_{1}(1)=0, \frac{dH_{1}}{d\xi}(-1)=1, \frac{dH_{1}}{d\xi}(1)=0$$
   
   

4. $H_{2}(\xi)$ associée à la valeur de la dérivée au second noeud ($\xi=+1$): i.e. telle que
   
   
   $$H_{2}(-1)=0, H_{2}(1)=0, \frac{dH_{2}}{d\xi}(-1)=0, \frac{dH_{2}}{d\xi}(1)=1$$
   
   

L'expression de ces polynômes est obtenue en utilisant les 4 conditions imposées. On obtiens ainsi les expressions suivantes:

$$N_{1}(\xi) = (1-\xi)^{2}(\frac{1}{2}+\frac{\xi}{4})$$

$$H_{1}(\xi) = \frac{1}{4}(1-\xi)^{2}(1+\xi)$$ (4.10)

$$N_{2}(\xi) = (1+\xi)^{2}(\frac{1}{2}-\frac{\xi}{4})$$

$$H_{2}(\xi) = -\frac{1}{4}(1+\xi)^{2}(1-\xi)$$

Ces fonctions de formes ont tracées sur la figure 5.2.

Figure 5.2: fonctions de forme hermitiennes

L'approximation $v^{h}$ sur l'élément de référence s'écrit ( en notant $h=x_{k}-x_{k-1}$):

$$v^{h}(\xi)=v_{k-1}N_{1}(\xi)+\frac{h}{2}(\frac{dv}{dx})_{k-1... ...(\xi)+v_{k}N_{2}(\xi)+\frac{h}{2}(\frac{dv}{dx})_{k}H_{2}(\xi)$$ (4.11)

On constate donc que les fonctions bases $\Phi_{k}$ sont égales aux fonctions de formes $N_{1}$ou $N_{2}$, mais que les fonctions de base $\Psi_{k}$ sont égales à $\frac{h}{2}H_{1}$ ou $\frac{h}{2}H_{2}$, car la dérivation dans l'élément de référence (i.e. par rapport à $\xi$) n'est pas égale à la dérivation physique (i.e. par rapport à x), et on a la relation:

$$\frac{dv^{h}}{dx}=\frac{dv^{h}}{d\xi}\frac{2}{h}$$

Sur l'élément $e_{k}$ les fonctions de base s'écrivent:

$$\Phi_{k-1}(x)=N_{1}(\xi), \Psi_{k-1}=\frac{h}{2}H_{1}(\xi), \Phi_{k}(x)=N_{2}(\xi), \Psi_{k}=\frac{h}{2}H_{2}(\xi)$$

4.2.2 Calcul du Lagrangien discret

Pour calculer le Lagrangien discret $\mathcal{L}(v^{h})$, on décompose les intégrales en une somme sur tous les éléments du maillage:

$$\mathcal{L}(v^{h})=-\frac{1}{2}\sum_{k=1}^{ne}\underbrace{\l... ..._{k}}qv^{h} dx\right)}_{W_{k}}+\underbrace{Fv^{h}(L)}_{W_{L}}$$

Ce Lagrangien contiens une somme de termes quadratiques $E_{k}$, de termes linéaires $W_{k}$ et un terme de bord $W_{L}$. Le terme $E_{k}$ correspond à l'énergie élastique de l'élément $e_{k}$, $W_{k}$ au travail des forces linéiques appliquées sur l'élément $e_{k}$ et le terme de bord $W_{L}$ au travail de la force appliquée en $x=L$.

D'un point de vue mécanique, nous avons décomposé le système mécanique (la poutre) en $ne$ morceaux, puis calculé l'énergie mécanique et le travail des forces externes sur chaque élément, pour ensuite en faire la somme.

4.2.2.1 énergie élastique élémentaire $E_{k}$

Sur un élément $e_{k}$, l'énergie élastique est calculée en effectuant le changement de variable (5.9) vers l'élément de référence et en utilisant l'expression (5.11) de $v^{h}$:

$$E_{k}=\frac{1}{2}\int_{-1}^{+1} EI_{zz}\left(\frac{4}{h^{2... ...\frac{dv}{dx})_{k}H_{2}(\xi)\right)\right)^{2} \frac{h}{2}d\xi$$

$E_{k}$ est donc une forme quadratique par rapport aux valeurs nodales $w^{k}=\left[v_{k-1}, (\frac{dv}{dx})_{k-1}, v_{k}, (\frac{dv}{dx})_{k}\right]^{t}$qui s'écrit:

$$\begin{eqnarray*} E_{k} & = & \frac{1}{2}\left(w^{k}\right)^{t} \mathbf{K}^{k}\... ...{dv}{dx})_{k-1}\\ v_{k}\\ (\frac{dv}{dx})_{k}\end{array}\right]\end{eqnarray*}$$

$\mathbf{K}_{k}$ est une matrice $(4,4)$ qui a pour expression:

$$\mathbf{K}^{k}=\frac{8EI_{zz}}{h^{3}}\left[\begin{array}{ccc... ...}}{4}\int_{-1}^{+1}\left(H_{2}\right)^{2}d\xi\end{array}\right]$$

C'est la matrice de raideur élémentaire. Cette matrice est symétrique, et on calcule ses coefficients en utilisant la définition des fonctions de forme (5.10). En utilisant les symétries de ces fonctions de forme, on ne calcule que 5 coefficients: $\mathbf{K}_{11}, \mathbf{K}_{12}, \mathbf{K}_{13}, \mathbf{K}_{22},\mathbf{K}_{24}$ et on déduit les autres par symétrie:

$$\mathbf{K}^{k}=\frac{EI_{zz}}{h^{3}}\left[\begin{array}{cccc... ...-6h & 12 & -6h\\ 6h & 2h^{2} & -6h & 4h^{2}\end{array}\right]$$ (4.12)

4.2.2.2 travail des forces élémentaire $W_{k}$

Sur un élément $e_{k}$, le travail des forces linéiques se calcule en effectuant le passage vers l'élément de référence:

$$W_{k}=\int_{-1}^{+1}q \left(v_{k-1}N_{1}(\xi)+\frac{h}{2}(\... ...+\frac{h}{2}(\frac{dv}{dx})_{k}H_{2}(\xi)\right)\frac{h}{2}d\xi$$

C'est une forme linéaire en fonction des valeurs nodales $w^{k}$, qui s'écrit:

$$\begin{eqnarray*} W_{k} & = & (w^{k})^{t}\mathbf{B}^{k}\\ & = & \sum_{i=1}^{4}w_{i}^{k} \mathbf{B}_{i}^{k}\end{eqnarray*}$$

$\mathbf{B}^{k}$ est un vecteur de dimension (4), qui a pour expression:

$$\mathbf{B}^{k}=\frac{h}{2}\left[\begin{array}{c} \int_{-1}^{... ...2}d\xi\\ \frac{h}{2}\int_{-1}^{+1}qH_{2}d\xi\end{array}\right]$$

C'est le vecteur des forces élémentaires, qui a pour expression dans le cas d'une force linéique $q$ constante:

$$\mathbf{B}^{k}=\frac{qh}{2}\left[\begin{array}{c} 1\\ +\frac{h}{6}\\ 1\\ -\frac{h}{6}\end{array}\right]$$ (4.13)

4.2.2.3 terme de bord $W_{L}$

le terme de bord se calcule facilement, puisque que la valeur de $v^{h}$ en $x=L$ vaut $v^{h}(L)=v_{nn}$. Il vient:

$$W_{L}=Fv_{nn}$$

4.2.3 Assemblage

Le Lagrangien discret s'écrit donc:

$$\begin{eqnarray*} \mathcal{L}(v^{h}) & = & -\frac{1}{2}\sum_{k=1}^{ne}E_{k}+\sum... ...thbf{K}^{k}w^{k}+\sum_{k=1}^{ne}(w^{k})^{t}\mathbf{B}^{k}+Fv_{nn}\end{eqnarray*}$$

En notant $\mathbf{X}$ le vecteur des $2nn$ valeurs nodales inconnues:

$$\mathbf{X}^{t}=\left[\begin{array}{cccccccc} v_{1} & (\frac{... ...dots & \ldots & v_{nn} & (\frac{dv}{dx})_{nn}\end{array}\right]$$

ce Lagrangien s'exprime sous la forme matricielle suivante:

$$\begin{eqnarray*} \mathcal{L}(\mathbf{X}) & = & -\frac{1}{2}\mathbf{X}^{t} \mat... ...{j=1}^{2nn}A_{ij}X_{i}X_{j}+\sum_{i=1}^{2nn}B_{i}X_{i}+FX_{2nn-1}\end{eqnarray*}$$

$\mathbb{A}$ est la matrice de rigidité globale du système obtenue par assemblage des matrices élémentaires $\mathbf{K}^{k}$. Elle est symétrique et a pour dimension $(2nn,2nn)$. $\mathbf{B}$ est le vecteur force global, obtenu par assemblage des vecteurs forces élémentaires $\mathbf{B}^{k}$ et a pour dimension $2nn$. Enfin le dernier terme du Lagrangien correspond à la contribution des conditions aux limites.

$\mathcal{L}(\mathbf{X})$ est une forme quadratique symétrique par rapport aux valeurs nodales $X_{i}$, et la condition de stationnarité de ce Lagrangien conduit aux $2nn$ équations de Lagrange:

$$\frac{\partial L(\mathbf{X})}{\partial X_{i}}=0$$

soit

$$-\sum_{j=1}^{2nn}A_{ij}X_{j}+B_{i}+\delta_{i,2nn-1}F=0\mbox{ pour }i=1,2nn$$

C'est un système linéaire de $2nn$équations:

$$\mathbb{A} \mathbf{X}=\mathbf{B}' \mbox{ avec }\mathbf{B}'=\mathbf{B}+\left[\delta_{i,2nn-1}\right]F$$

qu'il suffit de résoudre pour obtenir la solution approchée $v^{h}$. Ce système linéaire corresponds aussi à la formulation faible discrète de l'équation d'équilibre du système.

4.2.4 Applications

Dans le cas d'un maillage de $\Omega$ avec $ne=2$ éléments

$$\Omega=[0,\frac{L}{2}]\cup[\frac{L}{2},L]$$

la solution approchée possède donc au total $2nn=4$ degrés de liberté,correspondant aux valeurs nodales aux noeuds $x=\frac{L}{2}$ et $x=L$, les valeurs en $x=0$ étant imposées par les conditions aux limites.

$$v^{h}(x)=v_{1}\Phi_{1}(x)+(\frac{dv}{dx})_{1}\Psi_{1}(x)+v_{2}\Phi_{2}(x)+(\frac{dv}{dx})_{2}\Psi_{2}(x)$$

4.2.4.1 assemblage générique

De façon à effectuer l'assemblage de la matrice $\mathbb{A}$ et du second membre $B$ le plus simplement possible, nous allons tout d'abord considérer une solution approchée $v^{h}$ sans tenir compte des conditions aux limites, i.e. avec $6$degrés de liberté:

$$v^{h}(x)=v_{0}\Phi_{0}(x)+(\frac{dv}{dx})_{0}\Psi_{0}(x)+v_{... ..._{1}\Psi_{1}(x)+v_{2}\Phi_{2}(x)+(\frac{dv}{dx})_{2}\Psi_{2}(x)$$

Le vecteur inconnu $\mathbf{X}$ s'écrit

$$\mathbf{X}=\left[\begin{array}{cccccc} v_{0} & (\frac{dv}{dx... ...v}{dx})_{1} & v_{2} & (\frac{dv}{dx})_{2}\end{array}\right]^{t}$$

La matrice $\mathbb{A}$ est l'assemblage des 2 matrices élémentaires $\mathbf{K}^{1}$ et $\mathbf{K}^{2}$ associées aux 2 éléments du maillage:

$$\mathbb{A}=\left[\begin{array}{cccccc} K_{11}^{1} & K_{12}^{... ...1}^{2} & K_{42}^{2} & K_{43}^{2} & K_{44}^{2}\end{array}\right]$$

soit en utilisant l'expression (5.12) des matrices élémentaires avec $h=\frac{L}{2}$

$$A=\frac{EI_{zz}}{L^{3}}\left[\begin{array}{cccccc} 96 & 24L ... ... -24L\\ 0 & 0 & 24L & 4L^{2} & -24L & 8L^{2}\end{array}\right]$$

De même le second membre est l'assemblage des 2 seconds membres élémentaires $\mathbf{B}^{1}$ et $\mathbf{B}^{2}$:

$$\mathbf{B}=\left[\begin{array}{cccccc} B_{1}^{1} & B_{2}^{1}... ...{4}^{1}+B_{2}^{2} & B_{3}^{2} & B_{4}^{2}\end{array}\right]^{t}$$

ce qui donne en utilisant l'expression (5.13) du second membre élémentaire:

$$B=qL\left[\begin{array}{cccccc} \frac{1}{4} & \frac{L}{48} & \frac{1}{2} & 0 & \frac{1}{4} & \frac{L}{48}\end{array}\right]$$

4.2.4.2 Conditions aux limites

Les conditions aux limites du problème sont de 2 types:

1. les conditions de Dirichlet (conditions fortes): on remplace les équations 1 et 2 par
   
   
   $$v_{0}=0 \mbox{ et } (\frac{dv}{dx})_{0}=0$$
   
   

2. les conditions de Neuman (conditions faibles): on modifie le second membre pour ajouter le terme:
   
   
   $$\left[\delta_{i,2nn-1}\right]F$$
   
   

Le système linéaire s'écrit alors:

$$\frac{EI_{zz}}{L^{3}}\left[\begin{array}{cccccc} 1 & 0 & 0 &... ...t[\begin{array}{c} 0\\ 0\\ 0\\ 0\\ F\\ 0\end{array}\right]$$

La résolution de ce système linéaire fournit le vecteur solution approché:

$$X=[0, 0, \frac{1}{384}\frac{L^{3}\left(17qL+40F\right)}{EI... ...{EI_{zz}}, \frac{1}{6}\frac{L^{2}\left(qL+3F\right)}{EI_{zz}}]$$

La solution approchée s'écrit

$$\begin{eqnarray*} v^{h}(x) & = & \frac{1}{384}\frac{L^{3}\left(17qL+40F\right)}{... ...(x)+\frac{1}{6}\frac{L^{2}\left(qL+3F\right)}{EI_{zz}}\Psi_{2}(x)\end{eqnarray*}$$

c'est une fonction continue à dérivée continue et qui est un polynôme de degré 3 sur chaque élément.

La solution exacte de l'équation d'équilibre (5.1) associée aux conditions aux limites (5.2,5.3) est un polynôme de degré 4 qui s'écrit:

$$v_{ex}(x)=\frac{1}{24}\frac{qx^{4}}{EI_{zz}}-\frac{1}{6}\fra... ...}}{EI_{zz}}+\frac{1}{4}\frac{L\left(qL+2F\right)x^{2}}{EI_{zz}}$$

Figure: comparaison solution exacte et approchée pour $F=\frac{qL}{2}$ et $v_{ref}=\frac{qL^{4}}{EI_{zz}}$

On a tracé sur la figure (5.3) la solution exacte $v_{ex}$ et la solution approchée $v^{h}$ , ainsi que la norme de l'erreur $\vert v_{ex}-v^{h}\vert$ pour la cas $F=\frac{qL}{2}$. On constate que dans ce cas l'approximation par éléments finis du déplacement $v$ est excellente (les 2 courbes sont indiscernables) et on obtiens la valeur exacte aux noeuds du maillage. En particulier la méthode éléments finis prédit la déformation maximum en x=L:

$$v_{max}=\frac{L^{3}}{24EI_{zz}}(3qL+8F)$$

Figure 5.4: Comparaison entre la solution exacte et la solution par éléments finis pour l'angle de rotation $\theta$, le moment fléchissant $M_{z}$ et le cisaillement $V_{y}$

Sur la figure suivante (5.4), on a comparer les prédictions par éléments finis de l'angle de rotation, du moment fléchissant et de la contrainte de cisaillement avec la solution exacte. On constate une bonne prédiction de l'angle de rotation (dérivée première de $v$) et une assez bonne prédiction du moment fléchissant (dérivée seconde de $v$). Pour le cisaillement ( dérivée troisième de $v$), l'approximation par élément finis est constante par élément, alors que la solution exacte est linéaire. On ne prédit donc qu'une valeur moyenne par élément. On constate aussi que les conditions aux limites d'encastrement (ou conditions de Dirichlet) (5.2) sont vérifiées exactement par la solution approchée, alors que les conditions dynamiques (ou conditions de Neuman) ne sont vérifiées exactement pas:

$$EI_{zz}\frac{d^{2}v^{h}}{dx^{2}}(L)=-\frac{1}{48}qL^{2}\neq ... ...-EI_{zz}\frac{d^{3}v^{h}}{dx^{3}}(L)=F+\frac{qL}{4}\neq V_{y}=F$$

En particulier la valeur du cisaillement en $x=L$ correspond à la valeur moyenne du cisaillement dans le dernier élément puisque:

$$\frac{2}{L}\int_{\frac{L}{2}}^{L}V_{y}dx=F+\frac{qL}{4} \mbox{ car } V_{y}=F+q(L-x)$$