Un premier exemple

Table des matières

1. Un premier exemple#

Marc BUFFAT, dpt mécanique, Université Lyon 1

%matplotlib inline
from  numpy import *
import matplotlib.pyplot as plt

1.1. Etude d’une pale d’éolienne#

L’objectif de l’étude est de modéliser numériquement la déformation d’une pale d’éolienne sous l’effet de la rotation (force centrifuge). Pour cela nous utiliserons une approche mécanique, en construisant un modèle mécanique simplifié que l’on résoudra ensuite numériquement. Le modèle obtenu est équivalent à celui obtenu par une approche éléments finis générique, mais la démarche pour l’obtenir est spécifique au cas étudié.

eolienne

1.2. Démarche#

démarche de simulation numérique

  • modèle mathématique (EDP)

  • approximation (mécanique) de la solution

  • résolution numérique

  • analyse du résultat (erreur)

1.3. Modèle mathématique#

Le modèle utilisé est un modèle classique d’élasticité linéaire appliqué à une poutre de longueur \(R\) en traction sous l’effet de la rotation. Dans le repère lié à la pale on écrit

  • équilibre d’un élément de poutre de masse \(dm\) :

\[\Sigma \mbox{ forces de traction} = dm * \Gamma_{centrifuge} \]

  • \(\rightarrow\) équation d’équilibre (EDP) avec des conditions aux limites (CL) d’encastrement en \(r=0\) et de contrainte nulle en \(r=R\)

\[ \frac{\partial}{\partial r}\left(ES\frac{\partial u}{\partial r}\right)=\rho S\omega^{2}r\,\,\,\,\mbox{}{+C.L.\,}\,\,\, u(r=0)=0,\,\, ES\frac{\partial u}{\partial r}(r=R)=0 \]

1.4. Modèle discret d’une pale#

Au lieu de résoudre le problème continu précédent, on discrétise la poutre comme une série de petites poutres élémentaires, que l’on peut assimiler chacune à un système masse ressort, dont les propriétés mécaniques dépendent de la position (la masse et la raideur dépendent de la section).

modèle pale

1.5. Modèle discret d’un bout de poutre en traction#

On considère un élément de poutre modélisé par le système masse ressort suivant:

modèle mécanique

  • poutre = système masse + ressort

  • la raideur du ressort \(k\) est fonction du module Young \(E\):

\[ k=\frac{E\, S}{l}\]

  • bilan des forces appliquées:

    • force centrifuge sur les 2 masses

\[ \frac{m_1}{2} \omega^2 r_1 , \frac{m_2}{2} \omega^2 r_2\]

  • efforts extérieurs \(\vec{F}\) exercés par les 2 éléments de poutre voisins : \(f^1\) et \(f^2\)

  • équation du modèle sous forme matricielle pour chaque élément

\[\begin{split} \left[\begin{array}{cc} k & -k\\ -k & k \end{array}\right]\left[\begin{array}{c} u_{1}\\ u_{2} \end{array}\right]=\left[\begin{array}{c} f^1_{1} + \frac{m_1}{2} \omega^2 r_1\\ f^1_{2} + \frac{m_2}{2} \omega^2 r_2 \end{array}\right] \end{split}\]

  • regroupement (assemblage) des équations éléments par éléments

1.5.1. cas de 2 éléments#

Dans le cas de 2 éléments de masse et longueur identiques, on a 3 inconnues \(u_1,u_2,u_3\)

  • elts 1

\[\begin{split} \left[\begin{array}{cc} k & -k\\ -k & k \end{array}\right]\left[\begin{array}{c} u_{1}\\ u_{2} \end{array}\right]=\left[\begin{array}{c} f^1_{1} + \frac{m}{2} \omega^2 r_1\\ f^1_{2} + \frac{m}{2} \omega^2 r_2 \end{array}\right] \end{split}\]

  • elts 2

\[\begin{split} \left[\begin{array}{cc} k & -k\\ -k & k \end{array}\right]\left[\begin{array}{c} u_{2}\\ u_{3} \end{array}\right]=\left[\begin{array}{c} f^2_{2} + \frac{m}{2} \omega^2 r_2\\ f^2_{3} + \frac{m}{2} \omega^2 r_3 \end{array}\right] \end{split}\]

  • combinaison des équations pour éliminer les liaisons \(f^1_2 = - f^2_2\) $\( \left[\begin{array}{cc} k & -k & 0\\ -k & 2k & -k \\ 0 & k & -k \\ \end{array}\right]\left[\begin{array}{c} u_{1}\\ u_{2}\\ u_{3} \end{array}\right]=\left[\begin{array}{c} f^1_{1} + \frac{m}{2} \omega^2 r_1\\ m \omega^2 r_2 \\ f^2_{3} + \frac{m}{2} \omega^2 r_3 \end{array}\right] \)$

  • imposition des CL pour éliminer les efforts externes

  • encastrement en \(r=r_1\):
    l’équation \( u_1 = 0 \) remplace la 1ere équation

  • condition libre : contrainte nulle en \(r=r_3\)
    \(f^2_3 = 0\)

d’où le système linéaire 3x3 à résoudre

\[\begin{split} \left[\begin{array}{cc} 1 & 0 & 0\\ -k & 2k & -k \\ 0 & k & -k \\ \end{array}\right]\left[\begin{array}{c} u_{1}\\ u_{2}\\ u_{3} \end{array}\right]=\left[\begin{array}{c} 0\\ m \omega^2 r_2 \\ \frac{m}{2} \omega^2 r_3 \end{array}\right] \end{split}\]

généralisation avec N éléments

Par assemblage de N éléments on obtient un système linéaire de dimension \(N+1\), qui est équivalent à un modèle éléments finis!

1.5.2. modèle éléments finis discret#

Le système linéaire (modèle éléments finis) à résoudre est de la forme suivante:

\[\begin{split} \left[\begin{array}{ccccccc} {1} & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & {k}_{1}+{k}_{2} & -{k}_{2} & 0 & 0 & 0 & 0\\ 0 & -{k}_{2} & {k}_{2}+{k}_{3} & -{k}_{3} & 0 & 0 & 0\\ 0 & 0 & -k_{3} & k_{3}+k_{4} & -k_{4} & 0 & 0\\ 0 & 0 & 0 & -k_{4} & k_{4}+k_{5} & -k_{5} & 0\\ 0 & 0 & 0 & 0 & -k_{5} & k_{5}+k_{6} & -k_{6}\\ 0 & 0 & 0 & 0 & 0 & -k_{6} & {k}_{6} \end{array}\right]\left[\begin{array}{c} u_{1}\\ u_{2}\\ u_{3}\\ u_{4}\\ u_{5}\\ u_{6}\\ u_{7} \end{array}\right]=\omega^{2}\left[\begin{array}{c} 0\\ m_{2}r_{2}\\ m_{3}r_{3}\\ m_{4}r_{4}\\ m_{5}r_{5}\\ m_{6}r_{6}\\ m_{7}r_{7}/2 \end{array}\right] \end{split}\]

1.6. Application#

  • choix des paramètres

\(L=14 m\) , \(\omega=60 tr/min\) , \(M=800 kg\) , \(\rho = 8000 kg/m^{3}\) , \(E=200\,10^{9}\, N/m^{2}\) (acier)

  • raideur équivalente \(K \leftarrow ES/L\):

  • système linéaire équivalent $\( \mathcal{A} U = \mathcal{B} \)$

  • résolution à l’aide d’algorithme de type élimination de Gauss

  • étude de la précision en fonction du nombre d’éléments \(Ne\) : \(Ne=5,10,20\)

  • analyse: calcul de la contrainte en r=0

1.6.1. Algorithme#

  • définition des paramètres (variables)

  • calcul matrice \(\mathcal{A}\) et \(\mathcal{B}\) par assemblage (boucle sur les éléments)

  • résolution du système pour un nombre d’éléments \(N\) fixé

  • calcul de la déformée discrète et de la contrainte max

  • analyse de la précision en fonction de N

1.6.2. Programmation en Python#

  • définition des paramètres numériques

  • écriture de fonctions: assemblage, calcul de la solution (traction) en particulier pour l’étude paramétrique

""" modelisation d'une pale d'eolienne par elements finis (M. BUFFAT) """
import numpy as np
import matplotlib.pyplot as plt
# parametres
L=14.        # longueur
M=800.       # masse
rho=8000.    # densite
S=M/(rho*L)  # section equivalente
E=200e+09    # module d'Young
K=E*S/L      # raideur totale
w=60.        # rotation en tours/minute

def assemblage(Ae,ne):
    """ assemblage de la matrice A sur ne elts """
    A=np.zeros((ne+1,ne+1))
    for l in range(ne):
        A[l:l+2,l:l+2]=A[l:l+2,l:l+2]+Ae;
    # application de la C.L. en 0
    A[0,:]=0; A[:,0]=0; A[0,0]=1;
    return A

def traction(ne,w,L,M,K):
    """ calcul de la déformée en traction d'une pale en rotation w, de longueur L ,
        de masse M et de rigidité K  avec ne elts:
        renvoie le déplacement U à la position R et la tension en r=0.0  
    """
    omega=w*2*np.pi/60. # rotation en rd/s
    l=L/ne           # longueur d'un element
    m=M/ne           # masse elementaire
    k=K*ne           # raideur elementaire
    U=np.zeros((ne+1))    # deformation aux noeuds
    R=l*np.arange(0,ne+1) # position des noeuds
    F=omega**2*R          # force
    Ke=k*np.array([[1,-1],[-1,1]]) # matrice  elementaire de rigidite
    Me=m/2*np.array([[1,0],[0,1]]) # et de masse
    # calcul des matrices globales par assemblage
    Ka=assemblage(Ke,ne)
    Ma=assemblage(Me,ne)
    # resolution
    B=np.dot(Ma,F); B[0]=0; # calcul du second membre
    U=np.linalg.solve(Ka,B)
    # tension en O
    T0=np.dot(Ke[0,:],U[0:2])
    return U,R,T0

1.6.3. Calcul de la solution#

# etude de la solution en fonction du nbre d'elts
U1,R1,T1=traction(5,w,L,M,K)
U2,R2,T2=traction(10,w,L,M,K)
U3,R3,T3=traction(20,w,L,M,K)
# resultats
print("Tension  r=0: ",T1,T2,T3)
print("Deformee max: ",U1[-1],U2[-1],U3[-1])

Tension  r=0:  -221079.13858440143 -221079.13858440146 -221079.13858440102
Deformee max:  0.0014732713795264515 0.001451605623945181 0.0014461891850498607

1.6.4. Tracé et analyse du résultat#

# trace du résultat
plt.figure(figsize=(14,6))
plt.subplot(1,2,1)
plt.plot(R1,U1,'-x',lw=2,label="Ne=5")
plt.plot(R2,U2,'-o',lw=1,label="Ne=10")
plt.plot(R3,U3,'--',lw=2,label="Ne=20")
plt.legend(loc=0)
plt.title("Deformee pale ")
plt.subplot(1,2,2)
plt.plot(R1,U1-U3[::4],'-x',lw=1,label="Ne=5")
plt.plot(R2,U2-U3[::2],'-o',lw=1,label="Ne=10")
plt.legend(loc=0)
plt.title("Erreur ")
plt.draw()
plt.show()
../../_images/a991157019d7c00242dfa0aa38f2b8c228a2830c9cb4279937821e07e993e162.png

1.7. Résumé#

démarche appliquée

  1. Écriture du Modèle continu: EDP + CL: \(\mathcal{L}u=f\)

  2. Discrétisation du domaine en éléments (approche mécanique)

  3. Modélisation (mécanique) par éléments

  4. Construction du modèle discret:

  • système linéaire \(\mathcal{A}\{u^{h}\}=\mathcal{B}\)

  • assemblage de \(\mathcal{A}\) et \(\mathcal{B}\)

  • résolution

  1. Analyse du résultat: étude de la précision, calcul de la tension et de la déformée maximale

Contenu