Problème instationnaire (hyperbolique)

Contenu

8. Problème instationnaire (hyperbolique)#

8.1. Etude de la vibration d’une corde souple#

Considérons une corde tendue suivant l’axe \(Ox\), initialement au repos avec une longueur \(L\) et soumise à une tension \(T_{0}\)[1], que l’on déforme à l’instant \(t=0\). C’est le problème de la corde pincée d’un instrument de musique (clavecin), dont on se propose d’étudier la vibration (figure Fig. 8.1).

../_images/corde.png — Fig. 8.1 corde tendue souple au repos et à un instant t#

On suppose que la corde est sans raideur (ou souple), c’est à dire que la résultante des contraintes est la tension \(T(x)\) qui reste tangente à la corde. On néglige en particulier le moment de flexion qui apparaît si la corde possède une raideur (cas d’une poutre). On tiens cependant compte de l’élasticité de la corde, qui peut s’allonger proportionnellement à une variation de tension.

Considérons un élément de corde \(MM'\) de longueur \(ds\) (figure Fig. 8.1) situé en un point \(M\) de la corde. Au repos ce brin de corde \(MM'(0)\) est situé au point \(M(x,0)\) d’abscisse \(x\), d’ordonnée \(y=0\) et a pour longueur \(dx\) : \(M'=M'(x+dx,0)\)

A l’instant \(t\) il a subit un déplacement longitudinal \(u(x,t)\) et un déplacement transversal \(y(x,t)\). Le point \(M\) se trouve donc en \(x+u\) et \(y\): \(M=M(x+u,y)\) et le point \(M'\) en \(x+u+du\) et \(y+dy\): \(M'=M'(x+dx+u+du,y+dy)\). Le brin \(MM'\) fait un angle \(\theta(x,t)\) avec l’axe \(Ox\) et sa longueur \(ds\) vaut donc \(ds^{2}=(dx+du)^{2}+dy^{2}\), soit puisque les variables ne dépendent que de \(x\) et du temps \(t\):

\[ds=dx\sqrt{\left(1+\frac{\partial u}{\partial x}\right)^{2}+\left(\frac{\partial y}{\partial x}\right)^{2}}\]

L’angle \(\theta\) vérifie les relations suivantes:

\[\sin(\theta)=\frac{dy}{ds}=\frac{\frac{\partial y}{\partial x}}{\sqrt{\left(1+\frac{\partial u}{\partial x}\right)^{2}+\left(\frac{\partial y}{\partial x}\right)^{2}}}\,\,\mbox{ et }\,\,\cos(\theta)=\frac{dx+du}{ds}=\frac{1+\frac{\partial u}{\partial x}}{\sqrt{\left(1+\frac{\partial u}{\partial x}\right)^{2}+\left(\frac{\partial y}{\partial x}\right)^{2}}}\]

Pour de petits déplacements (\(du\ll dx,\,\,dy\ll dx)\), et à un instant \(t\) fixé on en déduit les relations suivantes:

\[\frac{\partial s}{\partial x}\simeq1+\frac{\partial u}{\partial x}+\frac{1}{2}\left(\frac{\partial y}{\partial x}\right)^{2}\]

\[\sin(\theta)\simeq\frac{\partial y}{\partial x}\,\,\mbox{ et }\cos(\theta)\simeq1\]

L’allongement relatif du brin \(MM'\) s’écrit:

\[\frac{ds-dx}{dx}=\frac{ds}{dx}-1\simeq\frac{\partial u}{\partial x}+\frac{1}{2}\left(\frac{\partial y}{\partial x}\right)^{2}\]

La tension dans la corde est donc la somme de la tension initiale \(T_{0}\) et d’une contrainte élastique, qui d’après la théorie de l’élasticité linéaire s’écrit:

\[T(x,t)=T_{0}+ES\left(\frac{ds}{dx}-1\right)\simeq T_{0}+ES\left(\frac{\partial u}{\partial x}+\frac{1}{2}\left(\frac{\partial y}{\partial x}\right)^{2}\right)\]

où \(E\) est le module d’Young de la corde et \(S\) sa section.

L’équation d’équilibre dynamique pour le brin \(MM'\) de masse volumique \(\rho\):

(8.1)#\[\begin{split}\begin{aligned} \rho Sdx\frac{\partial^{2}y}{\partial t^{2}} & = & \frac{\partial}{\partial x}\left(T\sin(\theta)\right)dx\simeq\frac{\partial}{\partial x}\left(T_{0}\frac{\partial y}{\partial x}+ES\left(\frac{\partial u}{\partial x}+\frac{1}{2}\left(\frac{\partial y}{\partial x}\right)^{2}\right)\frac{\partial y}{\partial x}\right)dx\\ \rho Sdx\frac{\partial^{2}u}{\partial t^{2}} & = & \frac{\partial}{\partial x}\left(T\cos(\theta)\right)dx\simeq\frac{\partial}{\partial x}\left(T_{0}+ES\left(\frac{\partial u}{\partial x}+\frac{1}{2}\left(\frac{\partial y}{\partial x}\right)^{2}\right)\right)dx\end{aligned}\end{split}\]

En ne conservant que les termes du premier ordre (petites oscillations), on obtiens un système de 2 équations découplées :

(8.2)#\[\begin{split}\begin{aligned} \rho S\frac{\partial^{2}y}{\partial t^{2}} & \simeq & T_{0}\frac{\partial^{2}y}{\partial x^{2}}\\ \rho S\frac{\partial^{2}u}{\partial t^{2}} & \simeq & \frac{\partial}{\partial x}(ES\frac{\partial u}{\partial x})\end{aligned}\end{split}\]

qui sont 2 équations des ondes, correspondant respectivement à des ondes de flexion ( \(y(x,t)\)) de célérité \(c_{y}=\sqrt{\frac{T_{0}}{\rho S}}\) et des ondes de compression (\(u(x,t)\)) de célérité \(c_{u}=\sqrt{\frac{ES}{\rho S}}\).

Dans la cas d’une corde de section \(S\) constante, ces 2 équations s’écrivent sous la forme d’une équation des ondes avec une célérité constante \(c_{0}\) (\(c_{0}=c_{y}\) ou \(c_{0}=c_{u}\)):

(8.3)#\[\frac{\partial^{2}v}{\partial t^{2}}=c_{0}^{2}\frac{\partial v^{2}}{\partial x^{2}}\]

Dans le cas d’une corde de section \(S(x)\) non constante (de section moyenne \(S_{0}\)), cette équation s’écrit:

(8.4)#\[\alpha\frac{\partial^{2}v}{\partial t^{2}}=c_{0}^{2}\frac{\partial}{\partial x}\left(\beta\frac{\partial v}{\partial x}\right)\]

avec \(\alpha=\frac{S}{S_{0}}\), \(\beta=1\) et \(c_{0}=\sqrt{\frac{T_{0}}{\rho S_{0}}}\) pour l’onde de flexion \(y(x,t)\) et \(\alpha=\frac{S}{S_{0}}\), \(\beta=\alpha\) et \(c_{0}=\sqrt{\frac{ES_{0}}{\rho S_{0}}}\) pour l’onde de compression \(u(x,t)\).

A cette équation on ajoutte des conditions aux limites de déplacement nul aux 2 extrémitées:

(8.5)#\[v(0,t)=0,\,\,v(L,t)=0\]

et des conditions initiales (CI) : déformation \(v_{0}(x)\) sans vitesse initiale

(8.6)#\[v(x,0)=v^{0}(x),\,\,\,\frac{\partial v}{\partial t}(x,0)=0\]

8.2. Formulation faible#

La formulation faible de l’équation des ondes (8.3) s’obtient classiquement en multipliant par une fonction test \(w(x)=\delta v\) (variation de \(v(x,t)\) à un instant t fixé), et en intégrant l’équation à \(t\) fixé sur le domaine de calcul:

\[\int_{0}^{L}\frac{\partial^{2}v}{\partial t^{2}}w\,dx=\int_{0}^{L}c_{0}^{2}\frac{\partial v^{2}}{\partial x^{2}}w\,dx\]

En intégrant par partie le membre de droite, et en utilisant les conditions aux limites de Dirichlet (CL) sur \(v\) (8.5), qui imposent des conditions homogènes sur la variation \(w\)

\[w(0)=0\,\,\,\mbox{ et }\,\,w(L)=0\]

on obtiens la formulation faible suivante:

(8.7)#\[\begin{split}\left\{ \begin{array}{c} \mbox{Trouvez }v(x,t)\,\mbox{ vérifiant les CI et CL et à chaque instant t}\\ \int_{0}^{L}\frac{\partial^{2}v}{\partial t^{2}}w\,dx+c_{0}^{2}\int_{0}^{L}\frac{\partial v}{\partial x}\frac{\partial w}{\partial x}\,dx=0\,\,\,\,\,\,\forall w(x)\,\mbox{ t.q. }\,w(0)=w(L)=0 \end{array}\right.\end{split}\]

Cette formulation faible correspond au théorème des travaux virtuels, qui traduit la condition de minimisation de l’action \(\mathcal{A}=\int_{0}^{\tau}\mathcal{L}(v,\dot{v})\,dt\), avec un Lagrangien:

\[\mathcal{L}(v,\dot{v})=\underbrace{\frac{1}{2}\int_{0}^{L}\left(\frac{\partial v}{\partial t}\right)^{2}dx}_{T=\mbox{énergie cinétique}}\,-\,\underbrace{\frac{1}{2}c_{0}^{2}\int_{0}^{L}\left(\frac{\partial v}{\partial x}\right)^{2}dx}_{U=\mbox{énergie potentielle}}\]

La formulation variationnelle associée s’écrit alors:

\[\begin{split}\left\{ \begin{array}{c} \mbox{Trouvez }v(x,t)\,\mbox{ vérifiant les CI et CL}\\ \mbox{minimisant l'action }\mathcal{A}=\int_{0}^{\tau}\mathcal{L}(v,\dot{v})\,dt \end{array}\right.\end{split}\]

8.3. Formulation éléments finis#

Pour construire une approximation \(v^{h}\) par éléments finis de la solution de (8.6), on utilise un maillage \(\mathcal{M}^{h}\) du domaine de calcul \(\Omega=[0,L]\) en \(ne\) éléments \(\left\{e_{k}=[x_{k-1},x_{k}]\right\}_{k=1,ne}\) associé à \(ne+1\) noeuds \(\left\{ x_{i}\right\} _{i=0,ne}\). On choisit ensuite une interpolation en espace de la solution. Pour notre étude nous utiliserons des éléments finis \(\mathcal{P}^{1}\). La solution approchée s’écrit alors comme une combinaison linéaire des fonctions de base \(\left\{ \Phi_{i}(x)\right\}\)associées à ce maillage \(\mathcal{M}^{h}\) et à l’interpolation \(\mathcal{P}^{1}\):

\[v^{h}(x,t)=\sum_{i=0}^{ne}v_{i}(t)\Phi_{i}(x)\]

En tenant compte des conditions aux limites (8.5), et des propriétés des fonctions de base \(\Phi_{i}(x_{j})=\delta_{ij}\):

\[v^{h}(0,t)=v_{0}(t)=0\,\,\mbox{ et }\,v^{h}(L,t)=v_{ne}(t)=0\]

l’approximation possède \(ne-1\) degrés de liberté et s’écrit:

\[v^{h}(x,t)=\sum_{i=1}^{ne-1}v_{i}(t)\Phi_{i}(x)\]

On note que, contrairement aux problèmes statiques étudiés précédemment, les coefficients \(v_{i}\) de la combinaison linéaire (valeurs nodales) dépendent maintenant du temps. Les fonctions tests \(w^{h}(x)\) étant des variations (à \(t\) fixé) de la solution \(v^{h}\), c.a.d. une combinaison des fonctions de base, on choisit comme fonctions tests les fonctions de base \(\Phi_{i}\). La formulation faible discrète s’écrit alors:

\[\sum_{j=1}^{ne}\frac{d^{2}}{dt^{2}}(v_{j})\int_{0}^{L}\Phi_{j}(x)\Phi_{i}(x)\,dx+c_{0}^{2}\sum_{j=1}^{ne-1}v_{j}(t)\int_{0}^{L}\frac{\partial\Phi_{j}}{\partial x}\,\frac{\partial\Phi_{i}}{\partial x}\,dx=0\,\,\,\forall i=1,ne-1\]

C’est un système de \(ne-1\) équations différentielles linéaires du second ordre qui s’écrit sous la forme matricielle:

\[\mathbf{M}\frac{d^{2}V}{dt^{2}}+c_{0}^{2}\mathbf{K}V=0\]

où \(\mathbf{M}\) est la matrice de masse, \(\mathbf{K}\) la matrice de rigidité et \(V\) le vecteur inconnu des valeurs nodales:

\[\mathbf{M}_{ij}=\int_{0}^{L}\Phi_{j}(x)\Phi_{i}(x)\,dx,\,\,\,\mathbf{K}_{ij}=\int_{0}^{L}\frac{\partial\Phi_{j}}{\partial x}\,\frac{\partial\Phi_{i}}{\partial x}\,dx,\,\,\,V_{i}=v_{i}(t)\]

Ce système s’écrit sous la forme suivante (noté le signe - devant la matrice \(\mathbf{A}\)):

\[\frac{dV^{2}}{dt^{2}}=-\mathbf{A}V\,\,\,\mbox{ avec }\mathbf{A}=c_{0}^{2}\mathbf{M}^{-1}\mathbf{K}\]

auquel on ajoute les conditions initiales (8.6). La matrice \(\mathbf{A}\) est une matrice symétrique et définie positive[2], puisque \(\mathbf{M}\) et \(\mathbf{K}\) sont des matrices symétriques définies positives. \(\mathbf{A}\)possède donc \(ne-1\) valeurs propres réelles positives, que l’on note \(\omega_{i}^{2}\).

La solution générale de ce système d’équations différentielles linéaires est donc la somme de \(ne-1\) solutions élémentaires où modes propres \(\Lambda_{i}\) :

\[V(t)=\sum_{i=1}^{ne-1}\alpha_{i}\Lambda_{i}e^{I\omega_{i}t}\]

les \(\Lambda_{i}\) étant les vecteurs propres associés aux valeurs propres \(\omega_{i}^{2}\) de la matrice \(\mathbf{A}\):

\[\mathbf{A}\Lambda_{i}=\omega_{i}^{2}\Lambda_{i}\]

D’un point de vue mécanique, la vibration du système est une combinaison de modes propres de vibration \(\Lambda_{i}\) associés à des fréquences propres \(\frac{2\pi}{\omega_{i}}\). On obtiens ainsi les premiers modes de vibrations de la corde.

8.4. Assemblage et calcul des modes propres#

Le calcul des matrices de masse et de raideur utilise la technique d’assemblage décrite dans les chapitres précédents.

On calcule tout d’abord une matrice de masse et de raideur élémentaire sur l’élément de référence, ce qui permet d’obtenir les matrices élémentaires de masse \(\mathbf{M}^{k}\) et de raideur \(\mathbf{K}^{k}\) sur un élément \(e_{k}\) de longueur \(h\):

\[\begin{split}\mathbf{K}^{k}=\frac{c_{0}^{2}}{h}\left[\begin{array}{cc} 1 & -1\\ -1 & 1 \end{array}\right],\,\,\,\mathbf{M}^{k}=h\left[\begin{array}{cc} \frac{1}{3} & \frac{1}{6}\\ \frac{1}{6} & \frac{1}{3} \end{array}\right]\end{split}\]

On effectue ensuite un assemblage élément par élément pour calculer les matrices globales de raideur \(\mathbf{K}\) et de masse \(\mathbf{M}\), puis on applique les conditions aux limites.

La prise en compte des conditions aux limites est effectuée en éliminant les lignes et les colonnes des 2 noeuds frontières (noeud 1 et noeud N).

Pour calculer les valeurs propres et vecteurs propres de la matrice \(\mathbf{A}=\mathbf{M}^{-1}\mathbf{K}\), on résout le problème équivalent de valeurs propres généralisés:

\[\mathbf{M}x=\lambda\mathbf{K}x\]

en utilisant une fonction de calcul de valeurs propres en spécifiant que les matrices sont symétriques définie positives.

8.5. Résultats#

Le calcul a été fait avec \(Ne=10\) éléments, soit \(9\) degrés de liberté. Le nombre de modes propres calculés est donc \(9\).

L’amplitude des modes propres de vibration calculés est tracé sur les figures ci-dessous pour les modes \(k=1\) \(k=6\) et \(k=8\) et comparé à la solution analytique \(\sin k\pi x\).

On constate que la précision diminue lorsque \(k\) augmente.

De même, on a tracé la courbe des fréquences propres calculés comparées aux fréquences analytiques \(k/2\)

../_images/frequence.png — Fig. 8.5 fréquences propres#

On constate la encore une diminution importante de la précision lorsque \(k\) augmente.

Conclusion: on vérifie le théorème de Shanon sur l’échantillonnage !