IN 103 : Introduction à Matlab Projet Vibrations de plaques

Vibrations des plaques minces

Le but de ce projet est de simuler les vibrations d’une plaque rectangulaire mince
lorsqu’elle est sujette à des efforts simples. Les vibrations de la plaque sont régies
par une équation aux dérivées partielles que l’on résolvera par projection modale.
L’excitation considérée sera un choc bref, que l’on approchera par un effort impul-
sionnel représenté par une distribution de Dirac. À l’aide des programmes, diverses
analyses seront menées afin de mettre en évidence le rôle des différents paramètres
physiques (analyse en fréquence, écoute du son produit par la plaque). Plusieurs simu-
lations seront faites, en variant les paramètres géométriques (longueur, largeur et épais-
seur), et en changeant la nature du matériau, ou uniquement son amortissement.

1 Définitions
On considère une plaque mince, dont l’épaisseur h est faible devant la longueur a et la largeur b, et
l’on supposera dans tout le projet que l’hypothèse des petites perturbations (H.P.P : petits déplacements
et petites déformations) est vérifiée. La plaque est située dans le plan (x, y), et le déplacement transverse,
qui est l’inconnue du problème, est noté w(x, y, t).

b y
w(x,y,t)

a h

Figure 1: Plaque rectangulaire, de dimensions a, b et d’épaisseur h.

L’équation qui gouverne les vibrations est donnée par le modèle de Kirchhof-Love:

D∆∆w + cẇ + mẅ = f (x, y, t). (1)

Dans cette équation, m est la masse surfacique (en kg.m−2) du matériau (m = ρh, où ρ est la masse
volumique), et D représente la raideur en flexion de la plaque:

Eh3
D= , (2)
12(1 − ν 2 )

1
où E est le module d’Young du matériau, et ν son coefficient de Poisson.
f (x, y, t) représente les efforts surfaciques extérieurs appliqués à la plaque (en N.m−2). Enfin, le
terme cẇ est un amortissement visqueux ad hoc que l’on a rajouté au modèle, et que l’on affinera par la
suite. Enfin, ∆ représente l’opérateur laplacien, de telle sorte que, en coordonnées cartésiennes, on a :
 4
∂4 ∂4

∂
∆∆ = + +2 2 2 (3)
∂x4 ∂y 4 ∂x ∂y
On considère des conditions aux limites de type “simplement supporté”:

∂2w
w=0 et =0 en x = 0 et x = a, (4a)
∂x2
∂2w
w=0 et =0 en y = 0 et y = b. (4b)
∂x2
Comme la plupart des problèmes de la physique, les mouvements sont régis par une équation aux
dérivées partielles (EDP) en temps et en espace, dont on ne connait pas de solution analytique. Des
méthodes d’approximation doivent donc être mises en œuvre pour la résolution. Dans le cadre de ce
projet, on utilisera la technique de projection modale que nous rappelons ci-après.
On s’intéresse dans un premier temps aux vibrations libres et non-amorties de la plaque, de telle sorte
que, dans l’équation (1), f (x, y, t) = 0 et c = 0. On recherche alors les mouvements harmoniques, i.e.
on cherche les solutions sous la forme : w(x, y, t) = Φ(x, y) exp(iωt), où ω est une pulsation donnée. Il
vient alors:
mω 2
∆∆Φ − Φ=0 (5)
D
On ajoute comme condition sur Φ(x, y) qu’elle doit vérifier les conditions aux limites (4). Un calcul
simple montre qu’une infinité de solutions vérifient cette équation. Elles sont indexées par deux nombres
entiers naturels, p et q, et s’écrivent:
 pπx   qπy 
Φp,q (x, y) = A sin sin . (6)
a b
A chacune de ces fonctions propres correspond une pulsation propre, qui est à la fréquence des oscilla-
tions du mode considéré:
r 
D p2 q 2

2
ωp,q = π + . (7)
m a2 b2
Une propriété essentielle des fonctions propres tient au fait qu’elles sont orthogonales vis à vis du produit
scalaire:
Z aZ b
< φ, ψ >= φ(x, y)ψ(x, y)dxdy. (8)
0 0

On peut alors choisir la constante A de l’équation (6), de telle sorte que la famille des fonctions {Φp,q }p,q=1...+∞
forme une base orthonormée par rapport au produit scalaire défini par (8). Soit finalement:
2  pπx   qπy 
Φp,q (x, y) = √ sin sin . (9)
ab a b
Une fois ce calcul réalisé, la propriété fondamentale d’orthonormalité est mise à profit pour projeter
la solution sur la base des modes propres, ce qui nous permet de découpler les variables temps et espace,

2
et de se ramener à un problème ne dépendant plus que du temps. On suppose que l’on peut décomposer
le déplacement inconnu w(x, y, t) selon:
+∞
X
w(x, y, t) = Xn (t)Φn (x, y), (10)
n=1

où n est un indice à deux entiers naturels permettant de numéroter les modes : n = (p, q), et Xn (t) est
appelé l’amplitude modale du mode considéré. En reportant l’expression (10) dans l’équation initiale
(1), puis en faisant le produit scalaire par un mode donné Φp (x, y), on obtient finalement un ensemble
d’équations aux dérivées ordinaires gouvernant la dynamique de chacune des amplitudes modales:
c
∀ p = 1... + ∞, Ẍp + Ẋp + ωp2 Xp = Fp (t) (11)
m
RaRb
où Fp (t) = m1 0 0 f (x, y, t)Φp (x, y)dxdy. Les avantages de cette formulation sont multiples. Tout
d’abord, le problème a été en quelque sorte résolu spatialement, puisque les seules inconnues restantes
désormais sont les amplitudes modales de chaque mode : Xp (t). De plus, les équations gouvernant la
dynamique des Xp (t) sont très simples : il s’agit d’équations d’oscillateurs linéaires amortis que l’on
sait résoudre dans de nombreux cas. Enfin, les oscillateurs sont découplés, si bien que l’on peut user de
la propriété fondamentale de superposition, propres aux problèmes linéaires : on peut calculer indépen-
damment la solution de chaque amplitude modale, puis recombiner toutes les solutions obtenues via (10)
afin de reconstruire le déplacement w(x, y, t).

Le but du projet est de programmer en Matlab les étapes rappelées en introduction, ainsi que de
résoudre les équations dynamiques (11) afin de pouvoir simuler la réponse de la plaque dans le cas de
forçages simples. On considèrera d’abord le cas d’un forçage impulsionnel, pour lequel on dispose d’une
solution analytique. La programmation permettra de faire ensuite un certain nombre d’analyse qui seront
exploitées : visualisation des oscillations de la plaque, écoute du son produit, analyse spectrale. Si le
temps le permet, on programmera ensuite le cas d’une excitation harmonique transitoire, pour laquelle
la résolution dynamique se fait par intégration temporelle directe des équations (11), ce qui correspond
aux cas généralement rencontrés en pratique, où l’on ne connait pas la solution analytique.
Une attention particulière sera portée au terme d’amortissement. Afin d’avoir des simulations plus
réalistes, le terme mc Ẋp sera différencié en fonction de la fréquence du mode. On écrira donc de manière
générique un amortissement différent pour chaque mode, et l’on adimensionera par la fréquence afin de
s’abstraire d’éventuels changements de dimensions géométriques de la plaque. Finalement, on remplac-
era donc systématiquement mc Ẋp par : 2ξp ωp Ẋp .

2 Calcul des modes propres

Le début du projet consiste à écrire quelques fonctions simples donnant les résultats fondamen-
taux dont on aura besoin pour l’analyse : modes propres et fonctions propres, selon les formules (9) et
(7). Il est demandé de mettre systématiquement (quand c’est nécessaire) comme arguments d’entrée les
paramètres physiques de la plaque : a, b, h, E, ν et m. Des valeurs numériques correspondants aux cas
d’une plaque en aluminium, en verre et en bois sont données section 6.

• Écrire un programme matlab qui prend comme arguments les paramètres physiques de la plaque,
et renvoit les déformées modales, ainsi que les pulsations propres. On pourra utiliser les fonctions
: mesh et meshgrid pour la représentation tridimensionelle, ainsi que sort pour trier les
modes par ordre croissant de fréquences.

3
 • Vérifier numériquement la propriété d’orthonormalité des modes propres.
On pourra utiliser la fonction : dblquad.

• Représenter les modes propres. Écrire une fonction prenant comme argument des valeurs données
d’amplitudes modales, et qui reconstruit le déplacement w(x, y) correspondant.
On pourra utiliser les fonctions : repmat et sum.

NB: A la fin de la séance 2, cette partie devra être achevée.

3 Réponse de la plaque à une excitation impulsionnelle

Le but de cette partie est de calculer la réponse vibratoire de la plaque à une excitation de type
impulsionnelle que l’on supposera, en première approximation, être une impulsion de Dirac. On suppose
de plus que cette impulsion est donnée en un point précis de la plaque noté (x0 , y0 ). La force extérieure
à considérer s’écrit donc:

f (x, y, t) = δ(x − x0 )δ(y − y0 )δ(t) (12)

Si bien que la dynamique à résoudre s’écrit:

∀ p = 1... + ∞, Ẍp + 2ξp ωp Ẋp + ωp2 Xp = Fp (t), (13)

Et Fp (t) vaut simplement:

1
Fp (t) = Φp (x0 , y0 )δ(t) = Fp δ(t). (14)
m
La solution de chaque équation modale est alors donnée simplement par:
Fp
Xp (t) = exp(−ξp ωp t)sin(ωpv t), (15)
ωpv

où la pulsation ωpv vaut:

q
ωpv = ωp 1 − ξp2 (16)

• Écrire une fonction prenant comme argument le point d’impact (x0 , y0 ) et renvoyant les valeurs
de la force à appliquer sur chaque oscillateur modal.

• Écrire une fonction calculant la réponse temporelle, selon la formule (15). On s’attachera à définir
proprement l’axe temporel.

• Représenter et animer la solution. Faites varier le nombre de modes pris en compte dans la simu-
lation.

4 Analyse des résultats

Afin d’analyser les résultats obtenus et vérifier leur pertinence physique, on pourra :

• Faire varier les dimensions géométriques (a, b, et h), afin de voir leurs influences sur les fréquences
propres.

4
 • Modifier les propriétés du matériau, en utilisant les valeurs données section 6.
• Tracer des spectrogrammes des signaux calculés (cf. fonction specgram).
• Ecouter le son produit par un point de la plaque. Pour ce faire, on utilisera le fait que, en première
approximation, la pression est proportionnelle à l’accélération de la plaque. On choisira donc un
point particulier de la plaque, pour lequel on a déjà calculé son déplacement. Par dérivation, on en
déduira l’accélération, que l’on écoutera grâce à la fonction soundsc.

5 Réponse de la plaque à une excitation sinusoïdale transitoire

On s’intéresse ici au calcul de la réponse temporelle de la plaque, lorsqu’on lui impose un forçage
sinusoïdal transitoire, appliqué en un point donné (x0 , y0 ). La force considérée s’écrit donc :
f (x, y, t) = δ(x − x0 )δ(y − y0 )s(t) sin(Ωt), (17)
où Ω est la pulsation extérieure imposée à la plaque, et s(t) une enveloppe choisie de telle sorte que :
– entre l’instant t = 0 et t = t1 , on augmente linéairement l’amplitude du forçage de 0 à une valeur A,
– entre t = t1 et t = tF , l’amplitude est gardée constante, égale à A.

t1 tF
Figure 2: Représentation de l’excitation s(t) sin(Ωt) considérée.

Afin de calculer la dynamique résultante de ce forçage, il est plus simple d’intégrer directement les
équations modales, qui auront la forme :
∀ p = 1... + ∞, Ẍp + 2ξp ωp Ẋp + ωp2 Xp = Fp s(t) sin(Ωt), (18)
où Fp = Φp (x0 , y0) / m.
• Programmer la résolution temporelle du système (18). Il est conseillé d’écrire plusieurs fonctions
pour simplifier la tâche : une fonction calculant le vecteur force en fonction du point d’impact, ainsi
qu’une autre renvoyant le membre de gauche de l’équation linéaire. Pour l’intégration temporelle,
on utilisera les fonctions : ode45 ou ode15s (cf. PC 6).
• Calculer les réponses pour différentes pulsations Ω et différents points d’impact sur la plaque, et
animer les solutions afin de mettre en évidence le régime transitoire. Commenter.
• Convergence modale : augmenter le nombre de modes pris en compte dans la simulation afin de
mettre en évidence les différences.

5
6 Valeurs numériques
6.1 Données géométriques
Il est fortement conseillé de laisser les paramètres géométriques a, b, et h comme données d’entrées
des programmes, afin de les laisser au choix de l’utilisateur. Ceci permettra de plus d’analyser les résul-
tats et d’interpréter les phénomènes physiques pour différentes valeurs. Cependant, pour commencer les
calculs, on pourra choisir par exemple:

a = 0.25 m
b = 0.35 m
h = 2 mm

6.2 Aluminium
On prendra les valeurs suivantes pour les constantes apparaissant dans le problème :

ρ = 2660 kg.m−3
E = 69 GPa
ν = 0.33

Pour l’amortissement, on pourra prendre par exemple: ∀ p, ξp = 1/ωp .

6.3 Verre
Pour le verre, on prendra les valeurs suivantes pour les constantes apparaissant dans le problème :

ρ = 2550 kg.m−3
E = 76 GPa
ν = 0.23

Pour l’amortissement, on pourra prendre par exemple: ∀ p, ξp = 0.0015.

6.4 Bois
On pourra prendre des valeurs qui se rapprochent du bois, cependant il faut bien avoir à l’esprit que
l’on a fait une hypothèse supplémentaire un peu forte, consistant à considérer que le bois est isotrope, ce
qui n’est pas vrai en réalité. Pour une modélisation plus précise du bois, il faudrait changer l’équation
de départ pour tenir compte de l’anisotropie.
A des fins d’écoute, on pourra cependant faire une simulation avec un matériau en “bois isotrope”,
avec les valeurs suivantes:
ρ = 388 kg.m−3
E = 77.3 GPa
ν = 0.26

Pour l’amortissement, on pourra prendre par exemple: ∀ p, ξp = 0.003.