Université de Pau et des Pays de l’Adour

TER

Impémentation des éléments finis pour

l’étude de l’effet de la géométrie sur
l’acoustique d’une pièce

Jean Dutheil
Tuteur :Sébastien Tordeux
M1 MMS

1er Février 2023 — 18 Mai 2023

Première partie
Introduction
La musique a toujours été un langage universel capable de transcender les
barrières culturelles et de créer des émotions intenses chez les auditeurs. Cepen-
dant, pour que cette expérience musicale soit véritablement immersive et capti-
vante, de nombreux éléments doivent être pris en compte, notamment l’acous-
tique de la salle de concert. La façon dont le son se propage et interagit avec
l’environnement peut grandement influencer la qualité de l’expérience auditive.
C’est là que la géométrie de l’espace entre en jeu. En effet, en optimisant l’ar-
chitecture des salles de concert, il est possible d’améliorer considérablement
l’acoustique et d’offrir ainsi une expérience sonore optimale aux spectateurs.
Dans cette optique, nous étudierons en 2 dimensions les conséquences sur le
son au sein de pièces de trois différentes géométries pour voir si la géométrie de
l’espace à un effet concret sur la propagation du son dans une pièce.

1
Deuxième partie
Étude analytique
1 Formulation forte
L’amplitude d’une onde sonore dans un domaine Ω est donnée par l’équation
différentielle suivante :

−∆u + k 2 u = f dans Ω (1)

Pour f une fonction source, k une constante valant le rapport entre la pul-
sation d’un son et sa célérité.
Avec les conditions aux bords de Fourrier-Robin

∂n u + u = α sur Γ (2)

Avec, α un nombre réel.

2 Formulation faible
On cherche à trouver une formulation faible de (5).
On a :

− ∆u + ku = f dans Ω
1 2
∀v ∈ H ,(∆u)v + k uv = f v dans Ω
Z Z Z
∀v ∈ H 1 , (∆u)v + k 2 uv = fv
Ω Ω Ω
Z Z Z Z
∀v ∈ H 1 , ∇.u∇.v − ∂n uv + k 2 uv = fv Formule de Green
ZΩ Γ
Z Z Ω Z Ω Z
∀v ∈ H 1 , ∇.u∇.v + k 2 uv + uv = f v + αv
Ω Ω Γ Ω Γ

Doncnotant :

∀u, v ∈ H 1 , a(u, v) = ∇.u∇.v + k 2

R R R R R
Ω Ω
uv + Γ
uv l(v) = Ω
fv + Γ
αv

On a la formulation faible du problème :

∀v ∈ H 1 , a(u, v) = l(v) (3)

2
3 Existence et unicité de la solution Faible
Pour montrer l’existence et l’unicité d’une solution à la formulation faible,
on va chercher à appliquer le théorème de Lax-Milgram.

1. H 1 est un espace de Hilbert

2. a est une forme bilinéaire, l une forme linéaire
3. lest linéaire : En effet, soit v ∈ H 1 (Ω)
Z
| f v| ≤ ||f ||L2 (Ω) ||v||L2 (Ω) ≤ ||f ||L2 (Ω) ||v||H 1 (Ω)
Ω
Z
| αv| ≤ Cα||v||H 1 (Ω) Continuité de la trâce
Γ

4. a est continue. En effet, soient u, v ∈ H 1

Z Z Z
2
|a(u, v)| = | ∇.u∇.v + k uv + uv|
ZΩ Ω
Z Γ
Z
2
≤| ∇.u∇.v| + k | uv| + | uv|
Ω Ω Γ
≤ ||∇u||L2 ||∇v||L2 + k 2 ||u||L2 ||v||L2 + C||u||H1 ||v||L2
Continuité de la Trace, inégalité de Holder
≤ (2 + k 2 + C)||u||H 1 ||v||H 1

5. a est coércive, en effet soit u ∈ H 1 :

Z Z Z
2 2
a(u, u) = (∇u) + k u + u2
2
Ω Ω Γ
≥ min(1, k 2 )||u||H 1

On peut donc appliquer le théorème de Lax-Milgram et déduire que la for-

mulation faible admet une unique solution dans H 1

Troisième partie
Éléments finies d’ordre k
4 Introduction
On considère un Hilbert V , muni d’un produit scalaire ⟨.⟩V et de la norme
associé ||.||V et on introduit le problème variationnelsuivant :

Trouver u ∈ V
∀v ∈ V, a(u, v) = l(v)

3
Avec a et l des formes respectivement bilinéaire et linéaire, continue et on a
coércive.
Le théorème de Lax-Milgram donne l’existence et l’unicité d’une solution a
ce problème variationnel. Le problème est que cet espace V est à priori de
dimension infinie est donc que ce problème n’est pas résoluble numériquement.

5 Dimension Finie
On cherche donc à approcher V par un sous-espace de dimension finie inclus
dans V qu’on notera Vh et on cherchera à résoudre avec les mêmes notations et
hypothèses que dansprécédente :

Trouver uh ∈ Vh
∀vh ∈ Vh , a(uh , vh ) = l(vh )

On appelle ce problème approché.

Theorème 1. Le problème approché admet une unique solution qui s’obtient

après résolution d’un système linéaire.
Démonstration. On sait que ce problème a une solution unique par le théorème
de Lax-Milgram , on note cette solution uh .
En effet, Vh est de dimension finie, donc complet, et avec la norme et le produit
scalaire induit par V , c’est un Hilbert.

Vh étant de dimension finie (supposons dim Vh = N ), on peut lui trouver une

base (ϕi )i∈[1,N ]
Ainsi :
N
X
∃Ui ∈ R, uh = Ui ϕi
i=1

On peut alorsécrire :
N
X
∀j ∈ [1, N ], a(ϕi , ϕj )Ui = l(ϕj )
i=1

Écrit sous formematricielle :

Ah U h = B h

Avec Ahi,j = a(ϕi , ϕj ), Ujh = Uj et Bjh = l(ϕj ) pour tout (i, j) ∈ [1, N ]2

4
Reste à montrer que Ah est définie positive :
Soit X ∈ RN :
N X
X N
⊤
X AX = xi a(ϕi , ϕj )xj
i=1 j=1
N
X N
X
= a( x i ϕi , x j ϕj )
i=1 j=1
N
X N
X
= a( x i ϕi , x i ϕi )
i=1 i=1
N
X
≥ α|| xi ϕi ||Vh a est coercive (de coefficient α)
i=1

PN
Le terme de gauche est nul si et seulement si la fonction i=1 xi ϕi est nulle et
comme les ϕi forment une base si et seulement si tout les xi sont nuls. On en
déduit que Ah est définie positive et que le système linéaire précédent admet
une unique solution.

6 L’espace P k de Lagrange
Pour,Ω un ouvert de R2 ,on propose ici une première possibilité de sous espace
de dimension finie de H 1 (Ω).

6.1 Triangulation
Dans l’ensemble du paragraphe, Ω désigne un ouvert de R2 .

La première étape consiste à ”découper” Ω pour obtenir un maillage trian-

gulaire ou triangulation conforme du domaine.

Notation 6.1. Pour Ω domaine de R2 , on note Πh (Ω) une triangulation de Ω.

Cette triangulation contient les triangles (ou éléments) Kh , les arrêtes Sh et
les sommets Nh , h représentant de pas du maillage.

Définition 6.1. On dit qu’un maillage Π est conformelorsque :

1. Les triangles sont non-dégénérés :les points d’un triangle ne sont pas
alignés
S
2. {Kh |Kh ∈ Πh } = Ω.
3. L’intersection de deux triangles distincs est soit vide, soit une arrête de
Sh soit un point de Nh .

5
Définition 6.2. On appelle triangle de base des éléments finies le triangle((0, 0), (0, 1), (1, 0)).
On le note Eb .
Définition 6.3. Pour chaque élément E ∈ Kh , on introduit TE la transforma-
tion affine suivante :
On note X1 = (x1 , y1 ), X2 = (x2 , y2 ), X3 = (x3 , y3 ) les trois nœuds de l’élément
E.
Et on définit T l’application qui d’un élément renvoie ces coordonnées barycen-
−−−−−→ −−−−−→
triques en X2 − X1 , X3 − X1 .
Propriétés 1.1. 1. Quelque soit l’élément E, celui-ci étant non-dégénéré,
TE est inversible.
2. Si x ∈ E alors TE (x) ∈ ((0, 0), (0, 1)(1, 0)).

6.2 Points d’interpolation

Le secondes étapes consiste à choisir des points d’interpolations sur les tri-
angles. Pour le faire, on définit des points d’interpolation sur le triangle de base
des éléments finis et on applique l’inverse de la transformation affine de chaque
élément pour trouver tous les points d’interpolations.

Définition 6.4. On définit ainsi :

Ike = {(m/k, n/k), m ∈ [0, k], n ∈ [0, k − m]}

On représente ici les points d’interpolations sur le triangle de base pour
k = 2, k = 3 ,k = 4

Figure 1 – Points ,d’interpolation pour k=2,3,4

6
7 Polynômes d’interpolation sur les éléments
7.1 Polynômes sur R2
Définition 7.1. On définit P l’ensemble des polynômes de R2 comme suit :

N X
M
2 X
P = {p : R2 → R, ∃M, N ∈ N, ∃(amn ) ∈ RN tels que :∀(x, y) ∈ R2 , p(x, y) = anm xn y m }
n=0 m=0

Notation 7.1. Dans la suite on introduit une numérotation des points d’in-
terpolations globale (I , J). Et pour chacun de ses indices, on note sI les coor-
données de ce point.
Définition 7.2. Et on définit Pk l’ensemble des polynômes de degrés au plus k
dans R2 comme suit :

2
P = {p : R2 → R, ∃M, N ∈ N, ∃(amn ) ∈ RN tels que :∀(x, y) ∈ R2 ,
N X
X M
p(x, y) = anm xn y m et(m + n > k) ⇒ amn = 0}
n=0 m=0

Définition 7.3. Enfin, pour un élément E, on définit Pk (E)

Pk (E) = {p : R2 → R, ∃q ∈ Pk , ∀(x, y) ∈ E, p(x, y) = q|E (x, y)}

7.2 Fonction Polynomial Par Element

On peut maintenant définir un espace de fonction définie tel que :

Vhk := {vh ∈ C 0 (Ω, ∈ Kh , vh|E ∈ Pk (E)} (4)

7.3 Les fonctions de formes

Soit I un point d’interpolations sur le triangle de base.
Alors, il existe (m, n) ∈ N2 tel que I = (m/k, n/k).
On pose alors :
1. φ1I (X) = Πm−1
i=0 (X − i/k)
2. φ2I (Y ) = Πn−1
i=0 (Y − i/k)
3. φ3I (Z) = Πk−n−m−1
i=0 (Z − i/k)
Enfin, on pose :

φI (X, Y ) = φ1I (X) × φ2I (Y ) × φ3I (1 − X − Y )

Propriétés 1.2. Soit I ∈ Ike , alors ∀J ∈ Ike différent de I, φI (J) = 0

7
Démonstration. En effet, on note I = (mI /k, nI /k) et J = (mJ /k, nJ /k).
Alors plusieurs cas sont possibles :
1. mJ < mI donc φ1I (J) = 0
2. nJ < nI donc φ2I (J) = 0
3. Dans le dernier cas, mJ + nJ > mI + nI donc k − mJ − nJ − 1 <
k − mI − nI − 1 et φ3I (J) = 0

Propriétés 1.3. Soit I ∈ Ike , alors les degrés de φI il vaut k.

Démonstration. Évidente il suffit d’ajouter les degrés des différents polynômes.

On peut voir ici un exemple de polynôme pour respectivement k = 2, 3.

Figure 2 – Points ,d’interpolation pour k=2,3,4

Ainsi, sur chaque élément on définit les fonctions de bases :

Définition 7.4. Pour un élément E ∈ Kh ,ie un point d’interpolation sur le
triangle élémentaire, la fonction de base sur E associé au point d’interpolation
ie est :
−1
∀x ∈ E : φE
ie (x) = φie ◦ TE (x) (5)
Enfin, on peut définir les fonctions de formes globales : Soit I un indice de
point d’interpolation global. On note KI l’ensemble des éléments contenant I et
pour chaque élément de KI , on note iE l’indice local du point d’interpolation I
dans E. Alors, on définit la fonction de forme associée à I :
X
ϕI = φE
iE
E∈KI

Donc pour tout uh ∈ Vh :

n
X
∃(ui ) ∈ R, uh = ui ϕi
i=0

8
8 Matrices
Pour I l’ensemble des points d’interpolation et pour i ∈ I, ϕi la fonction de
forme valant un en i, le problème discret devient, on note F l’ensemble des
fonctions de formes
Z Z Z Z Z
Trouver(ui )i∈I , tels que, ∀ϕj ∈ F , ui ∇ϕi ∇ϕj +ui ϕi ϕ j + ϕ i ϕj = f ϕj + αv
Ω Ω Γ Ω Γ

Le problème peut alors s’écrire comme solution d’un système linéaire dont
nous allons calculer les matrices

Quatrième partie
Implementation
9 Assemblage
L’algorithme d’assemblage est un algorithme qui crée la matrice de masse en
prenant en compte uniquement les contributions élémentaires de chaque fonction
de forme. En effet, il serait plus couteux de calculer l’intégrale sur l’ensemble
du domaine hors celle-ci est nulle pour beaucoup de couple de fonction et si elle
n’est pas nulle les fonctions sont nulles quasiment partout.
On peut décomposer chaque élément de la matrice A en fonction de l’apport
de chaque élément par linéarité de l’intégration :
X Z Z
A[i][j] = ∇ϕi ∇ϕj + ϕi ϕj
E∈K E E

On introduit, pour chaque élément E, aE la forme linéaire suivante :

Z Z
aE (ϕi , ϕj ) = ∇ϕi ∇ϕj + ϕi ϕj
E E

Ainsi :
XX
A[I][J] = a(ϕI , ϕJ )e⊤
I eJ
I J
XXX
A[I][J] = aE (ϕI , ϕJ )e⊤
I eJ
I J E
XXX
A[I][J] = aE (ϕI , ϕJ )e⊤
I eJ
E I J

On remarque que aE (ϕI , ϕJ ) est nul lorsque sI ou sJ ne sont pas des sommets
de E. Finalement, la somme sur tous les sommets du maillage se réduit à une
somme sur les points d’interpolation du triangle considéré !

9
 Il faut alors travailler localement sur tous les éléments, et il nous faut aussi
une fonction L2G qui, d’un élément est d’un indice de points d’interpolation du
triangle, renvoie l’indice global de ce point d’interpolation. Une fonction telle
que :
L2G(E, i) = I ⇐⇒ sI = sE i

Dans la suite donc on utilisera des lettres majuscules (I,J...) pour la numérotation
globale et des minuscules pour la numérotation locale.
Il vient alors :
XXX
E ⊤
A[I][J] = aE (φEi , φi )eL2G(E,i) eL2G(E,i)
E i j

On peut alors décrire l’algorithme d’assemblage :

Algorithm 1 Assemblage de matrices

1: Initialiser la matrice globale K avec des zéros
2: for chaque élément E do
3: for chaque nœud local i de l’élément e do
4: for chaque nœud local j de l’élément e do
5: calculer aE (φE E
i , φi )
6: l’ajouter à A[L2G(E, i),L2G(E, i)]

10 Matrice élémentaire
La matrice A peut être décomposée en deux matrices M , K matrice respec-
tivement de masse et de rigidité :

A = M + k2 K

Alors : Z
M [I][J] = ϕI ϕJ
Ω
Et : Z
K[I][J] = ∇ϕI ∇ϕI
Ω
Les contributions élémentaire peuvent aussi être décomposé de la même manière.
Z
ME [i][j] = φE E
i φj
E
Z
KE [i][j] = ∇φE E
i ∇φj
E

10
10.1 Calcul de la matrice élémentaire de masse
Soit E un élément fini,T E la transformation affine lui correspondant et
E
J sa jacobienne. Soient i, j deux indices locaux de points d’interpolations
dans cet élément. On note E e l’élément de base des éléments finis, le triangle
((0, 0)(0, 1)(1, 0)) : Alors :
Z
ME [i][j] = φEi φj
E
E
Z e
= | det T E | φei φej
E

Il suffit alors de calculer la matrice de masse sur l’élément de référence Me pour

toutes les avoir, en effet :
ME [i][j] = | det J E |Me [i][j]
Par exemple pour k = 1, la matrice élémentaire vaut :
 
2 1 1
1 
Me = 1 2 1
12
1 1 2

10.2 Calcul de la matrice élémentaire de rigidité

Avec les mêmes notations que sur le précédent paragraphe :
Z
KE [i][j] = ∇φE i .∇φj
E
E

On note en plus, ξ η les coordonnées barycentriques dans l’élément de base d’un

point. Or par la règle de la chaine, pour tout i :
" E# " # " ∂φe #
∂φi ∂ξ ∂η i
∂x ∂x ∂x ∂ξ
∂φE = ∂ξ ∂η ∂φei
i
∂y ∂x ∂y ∂η

On note BE la matrice de passage :

∇x,y φE e
i = BE ∇ξ,η φi

Inversement,
" ∂φe # " # " ∂φe #
i ∂x ∂y i
∂ξ
∂φei = ∂ξ
∂x
∂ξ
∂y
∂ξ
∂φei ∇ξ,η φei = JE ∇x,y φE
i
∂η ∂η ∂η ∂η

D’ou,
BE = (JE⊤ )−1
On a alors (parce que X.Y = X ⊤ Y :
Z
KE [i][j] = (φej )⊤ (BE
⊤
BE )φE
i
e

11
10.3 Second membre
10.3.1 Assemblage du second membre
Le second membre s’assemble comme le premier, mais plus facilement, voici
l’algorithme :

Algorithm 2 Assemblage de second membre

1: Initialiser le vecteur global K avec des zéros
2: for chaque élément E do
3: for chaque nœud local i de l’élément e do
4: calculer aE (f, φE
i )
5: l’ajouter à B[L2G(E, i)]

10.3.2 Second Membre élémentaire

De la même manière que précédemment, on définit :
Z Z e
BE [i] = E
f φi dxdy = |JE | f ◦ T −1 (ξ, η)φei (ξ, η)dξdη
E E

12
Cinquième partie
Code
11 Matrice Creuse
Voici une figure représentant les élément non nul d’une matrice : On re-

Figure 3 – Representation d’une matrice éléments finies

marque que celle-ci est vide en dehors d’une grande diagonale, il est donc inutile
d’enregistrer l’intégralité de la matrice. C’est même déconseillé, certaines de la
matrice crée sans ce format dépassait les 4 Gigas... Pour cela, on l’enregistre
au format creux grâce à la bibliothèque SciPy.sparse. Cette bibliothèque a aussi
l’avantage de proposer un solveur de matrice creuse efficace.

12 Géométrie
On veut étudier les effets de la géométrie d’une pièce sur la pression acous-
tique a l’intérieur. On veut comparer ces trois géométries : Le but est de com-
parer la pression acoustique à l’intérieur avec une onde libre, pour cela on se
munit d’une 4 géométrie, un disque plus grand que les pièces.
On pense sur ces géométries a ”forcer” des points de maillage pour pouvoir ef-
fectuer des mesures aux mêmes endroits sur les différentes structures.

13
 Figure 4 – Géométries, respectivement le rectangle, la salle 1, la salle 2

Figure 5 – Disque

13 Condition aux bords

On interprète les conditions aux bords de la manière suivante :

∂n u = ikZu

14
Avec i la base des nombres complexes, Z l’impédance acoustique du mur. Si Z=1
alors il n’y a pas de réflexion. Le système linéaire à résoudre est donc :

M + k 2 K − ikZMΓ = B

Ou MΓ et la matrice de masse au bord définit par :

Z
MΓ [I][J] = φI φJ
Γ

Elle se calcule facilement grâce à un algorithme d’assemblage non pas sur les
éléments mais sur les arrêtes du bord du domaine.

14 Source
On travaille avec la source suivante qu’on recentrera en un point si besoin :
5
1. S2(X, ϵ) = πϵ2 (1 − ||X||/ϵ) si X ∈ B(0, ϵ)
Suit une figure montrant deux sources sur le disque centrées en (−0.5, 0)(0.5, 0)

Figure 6 – Graphe de la source

15 Étude de convergence
Pour étudier la convergence, on travaille sur la géométrie rectangulaire et
on regarde si pour différents pas de maillages, la solution se rapproche de celle

15
 Figure 7 – Diagramme de convergence

avec le pas de maillage le plus petit.

16
Sixième partie
Resultats
16 Visualisation des simulations
Une première section dans laquelle on affiche les résultats. Elle a peu d’intérêt
qualitatif, mais elle permet plusieurs choses
1. Constater si le résultat est cohérent
2. Comprendre se qu’il se passe
On simule avec h le pas de maillage le plus petit que possible pour l’ordina-
teur, sans réflexions aux bords dans le disque et avec une impédance acoustique
de 2,6 pour les autres. On prend deux termes sources centrés en (1, 0.25)(2, 0.25).

Figure 8 – Graphes pour respectivement le rectangle, la salle 1 et la salle 2

17
 Enfin sur le disque :

Figure 9 – Graphes pour le disque

17 Amplitude a l’endroit de l’écoute.

On veut voir si la forme de la pièce influe sur la pression acoustique dans cer-
tains points. On effectue des mesures dans les suivants :(1, 1)(1.25, 1)(1.5, 1)(1.75, 1)(2, 2)
situer plus ou moins au centre et on prend deux sources en (1, 0.25)(2, 0.5). On
compare les résultats avec ceux obtenus dans le Disque sans réflexions aucune.

Figure 10 – Graphes pour respectivement le rectangle, la salle 1 et la salle 2

On remarque que les résultats sont plutôt similaires sans compter l’excrois-
sance à droite de ”Salle 1” je suppose que pas de maillage est responsable de
ceci..

18 Influence de l’impédance des murs

On effectue des simulations en faisant évoluer l’impédance des murs de la
pièce sur la pièce rectangulaire. On effectue des mesures au centre la pièce.

18
 Figure 11 – Graphes de l’influence de l’impédance

On remarque que l’impédance des murs à une forte influence sur l’intensité
reçue, en effet plus les murs sont ”réfléchissants” plus l’intensité est forte.

19 Réaction à différentes longueurs d’onde

Le but est de voir si les géométries renvoient une image correcte de la source
en domaine fréquentiel. On effectue pour chaque géométrie plusieurs simulations
pour F variant de 500 Hz à 15000Hz : Voici les graphes obtenus :

Figure 12 – Graphes pour respectivement le rectangle, la salle 1 et la salle 2

19
 Enfin, on peut comparer ces résultats avec ceux obtenus avec le disque.

Figure 13 – Graphes pour respectivement le rectangle, la salle 1 et la salle 2

On remarque que dans tous les cas, il semble y avoir un creux d’intensité vers
1020 Hz. Mais sinon, c’est la salle 1 qui semble au mieux reproduire les conditions
de propagation de l’extérieur, surtout dans les plus hautes fréquences.

20
Septième partie
Conclusions
20 Pertinence de la méthode
Globalement, la méthode des éléments finis semble pertinente pour la résolution
de ce genre de problème. Les résultats sont cohérents avec ce à quoi on pourrait
s’attendre.
Elle a cependant certains défauts :
1. 1) Les résultats sont dépendants du maillage, encore plus quand celui-ci
n’est pas très fin. Je pense que certaines difformités dans les résultats
viennent de la
2. 2) Bien que relativement efficace si on pousse le maillage à être un peu
petit, et avec un ordre de simulation grand, les assemblages de matrice
sont très longs. Dans notre cas, plus ou moins 2 h de calcul par géométrie.
Heureusement, elle a aussi de nettes avantages :
1. La résolution du système linéaire est très rapide ainsi que l’assemblage de
la matrice de second membre. Il est donc facile de changer les paramètres
de l’équation après assemblage des matrices et d’effectuer plusieurs simu-
lations avec les mêmes matrices.

21 Analyse des résulats

Il est difficile d’analyser les résultats, car ceci semble perturbé par le maillage
ou autre. Cependant, on peut conclure que la géométrie d’une pièce influe sur
l’intensité des ondes sonores à l’intérieur. De plus, la géométrie ”Salle 1” semble
avoir les meilleurs résultats fréquentiels, surtout dans les hautes fréquences. En-
fin, l’impédance des murs aussi affecte l’onde perçue. Les matériaux de construc-
tion sont donc primordiaux.

21