Modélisation et simulation de

systèmes discrétisés par la

méthode des éléments finis dans
le formalisme des Systèmes
Hamiltoniens à Ports : application
à la synthèse sonore

Raibaud Mathis
Master 2 ATIAM

Année universitaire : 2017-2018

Responsables du stage :
David Roze, Chargé de Recherche CNRS
Robert Piéchaud, Ingénieur de développement IRCAM
Résumé
Ce stage de fin d’études a pour but la modélisation et la simulation de systèmes de
dimension infinie non linéaires discrétisés par la méthode des éléments finis dans le
formalisme des Systèmes Hamiltoniens à Ports. Dans ce sujet, le modèle de corde non
linéaire de Kirchhoff-Carrier est étudié. La formulation faible du modèle non linéaire
permet d’en obtenir une version discrétisée approximant la solution par interpolation
polynômiale. En modifiant l’écriture du système, il est alors possible de décrire ce
problème dans le formalisme des SHP. A l’aide d’une méthode à passivité garantie
en temps discret, le bilan de puissance du problème non linéaire est assuré au cours
de la simulation numérique. Les résultats sont présentés et analysés pour différentes
valeurs d’excitations extérieures. Après une transcription du script dans un langage
compilé, les temps de calculs du modèle non linéaire en temps-réel sont discutés.
mots-clés : synthèse sonore, modélisation physique, systèmes non linéaires, méthodes
des éléments finis, Systèmes Hamiltoniens à Ports

Abstract
This master’s thesis aims at the modeling and the simulation of nonlinear infinite
dimensional systems discretized by the finite element method in the formalism of
Port Hamiltonian Systems. In this subject, the Kirchhoff-Carrier nonlinear string
model is studied. The weak formulation of the nonlinear model makes it possible to
obtain a discretized version approximating the solution by polynomial interpolation.
By modifying the writing of the system, it is then possible to adapt this problem in
the PHS’s formalism. Using a passivity method guaranteed in discrete time, power
balance of the nonlinear problem is ensured during the numerical simulation. The
results are presented and analyzed for different values of external excitations. After
a transcription of the algorithm in a compiled language, the computation times of
the nonlinear model in real time are discussed.
keywords : sound synthesis, physical modelling, nonlinear systems, finite element
method, Port Hamiltonian Systems

1
Remerciements
Je tiens à remercier mes encadrants David Roze et Robert Piéchaud pour avoir
accepté de me confier ce sujet, pour leur soutien et leurs conseils avisés tout au
long de ce stage. Je tiens également à remercier Thomas Hélie pour son accueil dans
l’équipe S3AM et pour sa grande implication auprès de l’ensemble des stagiaires.
Merci également à l’ensemble de l’équipe S3AM pour leur écoute et leurs avis durant
les présentations des travaux en cours en réunion d’équipe.
La formation ATIAM fut une expérience particulièrement riche, aussi bien du côté
professionnel que du côté humain. Pour ces deux raisons, j’aimerai remercier l’en-
semble des enseignants et des responsables de ce Master, avec une mention spéciale
pour Cyrielle Fiolet, dont le travail a été remarquable et sa disponibilité toujours
constante pour l’ensemble des élèves. Un grand merci évidemment à l’ensemble de
la promotion pour leur sympathie et leur convivialité.
Pour finir, je tiens à remercier sincèrement l’ensemble du personnel de l’IRCAM.
Ce haut lieu de la musique contemporaine et de la recherche scientifique est un
endroit exceptionnel et une réelle source de motivation. C’est une chance que d’avoir
pu travailler en ces murs tout au long de cette année.

2
Table des matières
1 Introduction 5
1.1 Contexte . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2 Etat de l’art . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.2.1 Musique et son : une histoire de mathématiques . . . . . . . . 6
1.2.2 Les différentes techniques de synthèse sonore numérique . . . . 6
1.2.3 Modélisation géométrique . . . . . . . . . . . . . . . . . . . . 9
1.2.4 Point de vue énergétique . . . . . . . . . . . . . . . . . . . . . 9
1.3 Notre approche . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

2 Équation des cordes 11

3 Méthode des éléments finis 13

3.1 Formulation faible . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
3.2 Discrétisation spatiale et interpolation polynômiale . . . . . . . . . . 13
3.3 Assemblage des matrices . . . . . . . . . . . . . . . . . . . . . . . . . 15
3.4 Conditions aux limites . . . . . . . . . . . . . . . . . . . . . . . . . . 16
3.4.1 Multiplicateurs de Lagrange . . . . . . . . . . . . . . . . . . . 16
3.4.2 Modification des éléments du système . . . . . . . . . . . . . . 17
3.4.3 Discussion sur les deux méthodes . . . . . . . . . . . . . . . . 18

4 Systèmes Hamiltoniens à Ports 19

4.1 Cas général . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
4.2 Structure de Dirac . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
4.3 Exemple : Oscillateur mécanique amorti . . . . . . . . . . . . . . . . 20

5 Modélisation d’une corde linéaire par la méthode des élément finis

dans le formalisme des Systèmes Hamiltoniens à Ports 22
5.1 Formulation faible du modèle linéaire . . . . . . . . . . . . . . . . . . 22
5.2 Approximation des énergies . . . . . . . . . . . . . . . . . . . . . . . 22
5.3 Ecriture dans le formalisme des SHP . . . . . . . . . . . . . . . . . . 23
5.4 Méthode numérique à passivité garantie . . . . . . . . . . . . . . . . 23
5.5 Simulation numérique . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

6 Extension au cas de la corde non linéaire 26

6.1 Modèle de Kirchhoff-Carrier . . . . . . . . . . . . . . . . . . . . . . . 26
6.2 Formulation faible du modèle non linéaire . . . . . . . . . . . . . . . 26
6.3 Approximation des énergies du modèle non linéaire . . . . . . . . . . 27
6.4 Gradient discret de l’hamiltonien non quadratique . . . . . . . . . . . 28
6.5 Résolution et simulation numérique . . . . . . . . . . . . . . . . . . . 29
6.6 Résultats . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

7 Vers un démonstrateur en temps-réel 32

7.1 Le framework JUCE . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
7.1.1 Application audio Standalone . . . . . . . . . . . . . . . . . . 32
7.2 Eigen : une bibliothèque d’algèbre linéaire . . . . . . . . . . . . . . . 33
7.3 Implémentations et résultats . . . . . . . . . . . . . . . . . . . . . . . 34

8 Conclusion et perspectives 35

3
9 Annexes 36
9.1 Annexe 1 : Force extérieure . . . . . . . . . . . . . . . . . . . . . . . 36
9.2 Annexe 2 : Comparaison des différentes bibliothèques d’algèbre li-
néaire en C++ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

Bibliographie 38

4
1 Introduction
1.1 Contexte
Ce stage, réalisé dans le cadre du Master 2 ATIAM (Acoustique, Traitement
du Signal et Informatique appliqués à la Musique, [1]), se déroule à l’Institut de
Recherche et Coordination Acoustique/Musique (IRCAM, [2]) dans l’équipe S3AM
(Systèmes et Signaux Sonores : Audio/Acoustique, instruMents, [3]) et est enca-
dré par David Roze, chargé de recherche CNRS, et Robert Piéchaud, ingénieur en
développement.
L’équipe S3AM est née de la fusion de deux équipes de l’IRCAM, anciennement
appelées "Acoustique Instrumentale" et "Systèmes et Signaux Sonores (S3)". De ce
mélange est née une équipe dont les travaux portent sur la modélisation, l’analyse,
le contrôle et la simulation de systèmes dynamiques tels que la voix, les instruments
de musique, les haut-parleurs ou encore les effets électroniques audio. Ces travaux
reposent sur les modèles physiques de systèmes dynamiques pour reproduire avec
réalisme le comportement des systèmes étudiés.
En 1988, l’IRCAM démarre le développement du logiciel Modalys [4] (ancienne-
ment nommé MOSAIC). Ce logiciel de synthèse sonore par modèles physiques per-
met de concevoir des instruments virtuels à partir de résonateurs linéaires (corde,
tube, plaque...) et d’excitateurs (archet, pincement, frappe...) qui interagissent en-
semble. Les objets sonores sont décrits de manière analytique ou par la méthode des
éléments finis dans un contexte modal. Actuellement, les travaux en cours envisagent
d’implémenter des résonateurs non linéaires dans Modalys.
Les recherches actuelles de l’équipe S3AM s’orientent depuis quelques années vers
l’utilisation du formalisme des Systèmes Hamiltoniens à Ports (SHP). Ce formalisme
permet l’interconnexion de systèmes multi-physiques (mécanique, électronique, chi-
mique, ...) tout en assurant la passivité du système global, même dans le cas non
linéaire. C’est dans ce contexte que l’intégration de résonateurs non linéaires dans
Modalys est envisagée. L’objectif de ce stage est de modéliser et simuler des systèmes
non linéaires approximés par la méthode des éléments finis dans le formalisme des
Systèmes Hamiltoniens à Ports. Le cas étudié dans ce projet est le modèle de corde
non linéaire de Kirchhoff-Carrier.
Ce document est organisé de la manière suivante. La partie 1.2 présente l’acous-
tique musicale d’un point de vue historique, décrit les différentes modélisations phy-
siques utilisées en synthèse sonore et discute de l’apport d’une vision énergétique
dans ces modélisations. La section 2 rappelle l’établissement de l’équation de la
dynamique d’une corde vibrante. La partie 3 détaille la méthode numérique des
éléments finis en l’appliquant au cas de la corde linéaire établie précédemment. La
partie 4 introduit le formalisme des Systèmes Hamiltoniens à Ports d’un point de vue
théorique et pratique. Des deux sections précédentes découle la modélisation d’une
corde linéaire discrétisée par éléments finis dans le formalisme des SHP en partie 5.
Cette modélisation est ensuite étendue dans la section 6 au modèle de corde non
linéaire de Kirchhoff-Carrier. Pour finir, l’implémentation d’un démonstrateur dans
un langage compilé est présentée et discutée.

5
1.2 Etat de l’art
1.2.1 Musique et son : une histoire de mathématiques

L’acoustique musicale, discipline mêlant sciences et musique, est née des travaux
du mathématicien Pythagore (VIème siècle avant J.-C.). Il découvre en effet qu’en
tendant une corde fixée à ses extrémités et en la bloquant en différents points spéci-
fiques avec son doigt, le pincement produits différentes notes agréables les unes avec
les autres. En mesurant ces longueurs, il constate que ces consonances apparaissent
lors de rapports entiers entre la longueur totale de la corde et la longueur tronquée
dû au positionnement de son doigt. Apparaissent ainsi l’octave, la quinte et la quarte
pour des rapports respectifs de 1/2, 2/3 et 1/3. Il faudra ensuite attendre le siècle
des Lumières et l’arrivée de Jean le Rond d’Alembert (1717-1783) [5] pour établir,
à partir des travaux d’Isaac Newton (1642-1727), de Jean Bernoulli (1667-1748)
et de Léonhard Euler (1707-1783), le mouvement d’une corde idéale en se basant
sur ses recherches sur les équations différentielles. La loi régissant la dynamique de
la corde, détaillée dans la partie 2, est alors appelée l’équation des ondes. Les re-
cherches théoriques du mathématicien Joseph Fourier (1768-1830) et expérimentales
de Hermann von Helmholtz (1821-1894) sur l’analyse des sons complexes ont ensuite
montré qu’à tout son musical de hauteur donnée est associé un timbre qui résulte
de la superposition au son fondamental d’une série d’harmoniques. L’ensemble de
ces découvertes montre ainsi le lien étroit entre ces deux domaines que sont la mu-
sique et la science. Les inventeurs du XXème siècle prennent alors parti de ce lien
dans leurs innovations et proposent de nouveaux instruments capables d’imiter les
instruments traditionnels. On peut citer entre autres le Thérémine inventé en 1919
par le Russe Lev Sergueïevitch Termen (connu sous le nom de « Léon Theremin »),
les Ondes Martenot par le Français Maurice Martenot présentées au public en 1928,
l’Ondioline par le Français Georges Jenny en 1941, le Moog Modular par Robert
Moog en 1964 et bien d’autres.

1.2.2 Les différentes techniques de synthèse sonore numérique

La seconde moitié du XXème sicèle, soit deux cent ans après la découverte de
l’équation des ondes par d’Alembert, les ordinateurs se démocratisent et permettent
ainsi à des scientifiques de réaliser de nombreux calculs en un temps record. Grâce
à cette nouvelle technologie, les premiers essais de synthèse sonore par ordinateur
voient le jour. Nous différencions ici les synthèses dites par approche "signal" de
celles par approche "physique". La première catégorie analyse les signaux obtenus
par une source et cherche à comprendre comment les reproduire. La deuxième tente
de comprendre et de simuler le comportement physique des sources étudiées. Les
exemples de synthèse par approche "signal" les plus connus sont : synthèse additi-
ve/soustractive, synthèse AM/FM, synthèse par tables d’ondes, synthèse granulaire.
Les approches "physiques" les plus employées en synthèse sonore sont détaillées ci-
dessous.

Synthèse par assemblage d’éléments physiques élémentaires

Le premier paradigme de modélisation physique est apparu à la fin des années

1970 et au début des années 1980. Il fut proposé par Cadoz [6] en collaboration avec
Luciani et Florens [7] et est à l’origine des environnements CORDIS et CORDIS-
ANIMA [8]. Il se base sur un ensemble de masse-ressort connectés les uns aux autres

6
et qui peuvent interagir grâce à des forces de contact non-linéaires. Cette représenta-
tion permet alors à l’utilisateur de créer des réseaux tel qu’il le souhaite. Cependant,
il est difficile de modéliser un système en particulier (corde, poutre, plaque) puisque
le réseau ne prend en compte ni la géométrie ni les caractéristiques physiques du
problème.

Synthèse modale

La synthèse par approche modale, développée initialement par Adrien [9], se base
sur le fait que le mouvement d’un objet physique peut être décomposé en une somme
infinie de sinusoïdes amorties multipliées par un facteur nommé déformée modale.
Afin de simuler numériquement le déplacement en un point de l’objet, le nombre de
modes est tronqué à N. Soit, pour un cas à une dimension,
N
X
w(x, t) = w̃(t)ψn (x), (1.1)
n=0
N
X
w(x, t) = an (t)e(−dn t) sin(2πfn t)ψn (x) (1.2)
n=0

avec w(x, t) le déplacement au point x et au temps t, et où ψn (x), an (t), dn et

fn correspondent respectivement à la déformée modale, à l’amplitude, à l’amortis-
sement et à la fréquence du mode n. ψn , dn et fn dépendent des caractéristiques
physiques de l’objet (géométrie, matériau) et de ses conditions aux limites. Leurs
expressions peuvent être obtenues de manière analytique ou par la méthode des élé-
ments finis lors de cas à géométrie complexe en résolvant un problème aux valeurs
propres généralisé. L’amplitude temporelle an (t) dépend quant à elle des conditions
initiales du problème et/ou de l’excitation extérieure. Le logiciel Modalys, développé
à l’IRCAM, s’appuie sur ce formalisme et permet en plus des interactions entre les
différents objets modélisés afin de concevoir des instruments virtuels utilisables en
temps réel.

Guides d’ondes numériques

De manière générale, la modélisation physique est réaliste mais nécessite de nom-

breux calculs pour s’approcher de la réalité. Au cours des années 1980, la synthèse
par guides d’ondes numériques, inspirée de l’article MSW [10] et développée par
Julius Smith [11], s’est imposée dans le domaine de la synthèse sonore par modèle
physique. Cette nouvelle synthèse, à l’inverse, a réussi à allier réalisme et faible coût
de calcul, particulièrement dans le cas de modélisation 1D (corde, tube). L’idée est
la suivante : le mouvement du milieu peut être décomposé en une somme de deux
ondes se propageant dans des directions opposées. Si c est la vitesse de l’onde se
propageant dans le milieu, la solution de l’équation peut être discrétisée en temps
avec une période T et avec un intervalle X en espace tel que X = cT . Les positions
et temps discrets sont notés x = mX et t = nT . Le déplacement s’exprime alors :

w(m, n) = w+ (m − n) + w− (m + n), (1.3)

où w+ et w− correspondent respectivement à l’onde se propageant vers la droite et
vers la gauche au point m et au temps n. Les guides d’ondes numériques repose

7
ensuite sur des lignes à retard pour prendre en compte la géométrie du système et
c’est à ce niveau que le coût de calcul est considérablement diminué. En effet, à
partir de la longueur parcourue par l’onde et sa vitesse, il est possible de convertir le
temps de propagation de cet aller-retour en nombre d’échantillons, ce qui ne nécessite
aucune opération coûteuse en temps de calcul pour l’ordinateur. Pour prendre en
compte les autres influences physiques sur le son (amortissement, dispersion), la
synthèse par guide d’ondes numérique s’appuie sur des techniques de traitement
du signal et particulièrement de filtrage. Par exemple, les réflexions aux limites du
système peuvent être modélisées par un filtre passe-bas, la dissipation par un gain
|g| < 1, les pertes dépendantes des fréquences par des filtres à phase linéaire.

Synthèse par la méthode des différences finies

La méthodes des différences finies est une méthode numérique permettant d’ap-
procher les solutions d’équations aux dérivées partielles. Elle se base sur une dis-
crétisation de l’espace et/ou du temps en un maillage régulier ou non. Le but est
d’approximer les valeurs des fonctions inconnues en certains points à partir des va-
leurs des points voisins. Pour ce faire, les approximations des dérivées premières et
secondes sont obtenues à partir de la formule de Taylor-Young au premier et second
ordre. Soit u(t) une fonction dépendant du temps t et δt un pas de temps, nous don-
nons ici les expressions des dérivées en fonction du type de différenciation choisie.
Pour la dérivée première,

u(t + δt) − u(t)

= u0 (t) + O(δt), Différenciation décentrée avancée (1.4)
δt

u(t) − u(t − δt)

= u0 (t) + O(δt), Différenciation décentrée retardée (1.5)
δt

u(t + δt) − u(t − δt)

= u0 (t) + O(δt), Différenciation centrée (1.6)
2δt
Pour la dérivée seconde,
u(t + δt) − 2u(t) + u(t − δt)
= u00 (t). (1.7)
δt2
Afin d’illustrer la simplicité de cette méthode, l’exemple de l’oscillateur harmo-
q
k
nique est étudié. Dans le cas où l’oscillateur possède une pulsation propre w0 = m ,
où k et m correspondent respectivement à la raideur du ressort et à la masse, un
raisonnement mécanique permet d’établir l’équation :

ẅ(t) + ω0 w(t) = 0. (1.8)

En approximant l’accélération ẅ(t) par l’équation (1.7),
w(t + δt) − 2w(t) + w(t − δt)
+ ω0 w(t) = 0. (1.9)
δt2
En mettant tous les termes dépendant de t et (t − δt) à droite, un schéma explicite
apparaît :

8
 w(t + δt) = w(t) 2 − ω0 δt2 − w(t − δt).
 
(1.10)
L’équation différentielle peut finalement être simulée numériquement, avec tn = nt
et δt = tn+1 −tn = n+1−n = 1. Les conditions initiales sont utilisées pour initialiser
les premières valeurs de w[n] et w[n − 1]. Finalement, à temps discret,

w[n + 1] = w[n] 2 − ω0 δt2 − w[n − 1].

(1.11)
Stefan Bilbao a étendu cette technique à de nombreux cas d’acoustique musicale et
les détaille dans [12] avec une analyse des conditions de convergence.

1.2.3 Modélisation géométrique

Dans le domaine de l’industrie mécanique, afin de prédire la réponse d’une struc-

ture à des efforts extérieurs, la modélisation d’objets à géométrie parfois complexe
est primordiale. Pour ce faire, la structure est approximée numériquement par la
méthode des éléments finis (cf. section 3). Cette méthode permet de décomposer
l’objet en un ensemble fini de sous-domaines délimités par des nœuds. Dans le cas
d’un système dépendant de l’espace et du temps, l’approximation par la méthode
des éléments finis conduit au système suivant :

MẄ + KW = F. (1.12)
avec M, K les matrices de masse et de raideur de dimension n × n, F le vecteur de
force de dimension n et W le vecteur de déplacements des nœuds de dimension n, où
n est le nombre de nœuds. L’approximation en espace est alors effectuée mais pas
celle en temps. Différentes techniques de discrétisation sont alors employées pour
résoudre le problème temporel [13, § 9] telles que la méthode des différences finies,
la technique de Houbolt [13, § 9.2.2], la méthode de Newmark [13, § 9.2.3] ou encore
la méthode de Bathe [13, § 9.2.4]. Cependant, la stabilité de ces méthodes n’est pas
garantie au cours du temps et dépend de différents paramètres propres à chaque
méthode.

1.2.4 Point de vue énergétique

Au XVIIème siècle, Isaac Newton exposait sa théorie de la mécanique classique.

La seconde loi exprime le fait que, dans un référentiel galiléen, la somme des forces
appliquées à un objet ponctuel est égale au produit de la masse de l’objet par
son vecteur accélération. En 1833, William Rowan Hamilton introduit la mécanique
Hamiltonienne, à partir des travaux de Joseph Lagrange. Ce formalisme est une im-
portante reformulation de la mécanique classique de Newton basé sur une approche
énergétique, menant plus tard aux formulations de la physique statistique et de la
mécanique quantique. L’avantage de ce formalisme est de pouvoir évaluer la quantité
d’énergie du système.
L’emploi de la mécanique hamiltonienne dans la synthèse sonore s’est avéré être
un atout majeur. En effet, en prenant en compte la conservation de l’énergie dans
les différentes techniques de synthèse sonore étudiées précédemment, il est possible
d’assurer la stabilité des simulations et ce, même dans le cas de systèmes non li-
néaires [14, 15, 16]. A partir de cette formulation de la mécanique, il est possible
de connecter plusieurs sous-systèmes pour en créer un global tout en assurant la
conservation de l’énergie aux points d’interconnexion. Le formalisme des Systèmes

9
Hamiltoniens à Ports [17, 18], introduit dans les années 1990 par Arjan Van der
Schaft et Bernhard Maschke comme une extension de la mécanique hamiltonienne,
repose sur cette idée. Employé en premier lieu dans le domaine de l’automatique
puis plus tard dans de nombreux domaines tels que la mécanique, l’électronique ou
encore la chimie, ce formalisme a permis de résoudre des problèmes de dimension
finie linéaire ou non linéaire pouvant inclure des phénomènes de dissipation (ex : os-
cillateur harmonique amorti, circuit RLC). Pour résoudre des systèmes de dimension
infinie (comprenant des équations aux dérivées partielles), une réduction d’ordre du
système est nécessaire. Celle-ci peut être assurée par différentes techniques, comme
par exemple la méthode des éléments finis mixtes [19], la méthode pseudo-spectrale
[20], les fonctions de Casimir [21], le formalisme modal [22] ou plus récemment la
méthode des éléments finis [23]. Bien qu’ayant des applications plutôt orientées vers
l’industrie, le formalisme des SHP a été utilisé lors de deux travaux de thèse [24, 25]
au sein de l’équipe S3AM de l’IRCAM à des fins musicales et continue d’être employé
dans le cadre de différents sujets de thèse et de stage (simulation d’ondes Martenot,
reproduction du conduit vocal, correction de haut-parleurs...).

1.3 Notre approche

Nous avons vu dans la partie 1.2.2 que le logiciel Modalys est basé sur le for-
malisme modal. Chaque objet créé possède ses fréquences propres et ses déformées
modales. Dans le cas de géométrie simple (corde homogène fixée à ses extrémités,
plaque mince simplement supportée sur ses bords, ...), ces variables sont calculées de
manière analytique. Lorsque l’objet conçu présente une géométrie complexe, la mé-
thode des éléments finis est employée pour les calculer. Cependant, les résonateurs
(corde, plaque, tube) implémentés dans Modalys sont à l’heure actuels linéaires. Le
formalisme modal n’est pas adapté à la modélisation de résonateur non linéaires. Le
passage par d’autres formalismes tels que les séries de Volterra est nécessaire mais
il n’est pas possible de déterminer a priori la stabilité du système. D’un autre côté,
c’est exactement le point fort du formalisme des Systèmes Hamiltoniens à Ports.
Notre approche consiste alors à prendre avantage de cette caractéristique pour pou-
voir implémenter des résonateurs non linéaires dont la stabilité de la simulation est
d’avance garantie. La technique de discrétisation spatiale retenue est la méthode
des éléments finis. Plusieurs raisons à ce choix : la première est que la méthode des
éléments finis nous permet de considérer des caractéristiques physiques non homo-
gènes (comme par exemple une corde de section variable, ou de matériau changeant
en fonction de la longueur) et la seconde est que cette méthode repose sur une dé-
composition du domaine en un nombre fini d’éléments et se rapproche ainsi plus de
la notion de sous-objets interconnectés (comme dans le formalisme des SHP). Afin
d’illustrer cette nouvelle méthode, le système considéré est une corde fixée à ses ex-
trémités et excitée par une force extérieure. Nous étudierons dans un premier temps
le cas linéaire et étendrons la méthode au modèle non linéaire de Kirchhoff-Carrier.

10
2 Équation des cordes
Cette partie a pour but de rappeler l’équation des cordes obtenue par un raison-
→
−
nement mécanique. La corde, de masse linéique µ et tendue par une tension T , est
représentée sur la figure 1.

Figure 1 – Bilan des forces appliquées à un tronçon de corde, source [26]

A partir du Principe Fondamental de la Dynamique, le mouvement d’ une corde

dans un plan (figure 1) est régi par les équations :

∂ 2w ∂ 2ξ
ρA = T sin(θ), ρA = T cos(θ), (2.1)
∂t2 ∂t2
∂ξ
avec ρ la masse volumique, A l’aire de la section, θ = arctan[ ∂w
∂x
/(1+ ∂x )], T la tension
totale, et où w et ξ correspondent respectivement aux mouvements transverse et
longitudinal de la corde. En reportant l’expression de θ dans les équations (2.1) et
en posant µ = ρA la masse linéique de la corde,
" #
∂ 2w ∂ T (∂w/∂x)
µ 2 = p , (2.2)
∂t ∂x (∂w/∂x)2 + (1 + (∂ξ/∂x))2
" #
∂ 2ξ ∂ T (1 + (∂ξ/∂x))
µ 2 = p , (2.3)
∂t ∂x (∂w/∂x)2 + (1 + (∂ξ/∂x))2
p

où T − T0 = EA (∂w/∂x)2 + (1 + (∂ξ/∂x))2 − 1 , avec T0 la tension statique
et E le module d’Young de la corde. Dans ce projet, les ondes longitudinales sont
négligées, soit ∂ξ/∂x ≈ 0. La dynamique de la corde s’exprime finalement [27] :
" p
#
∂ 2w ∂ T0 + EA( 1 + (∂w/∂x)2 − 1) (∂w/∂x)
µ 2 = p . (2.4)
∂t ∂x 1 + (∂w/∂x)2

11
Dans le cas où les coefficients d’amortissement fluide et structurel, respectivement α
et κ, sont introduits et où la corde est excitée par une force extérieure f , l’équation
des cordes (2.4) devient :

" p
#
∂ 2w ∂w ∂ 3w ∂ T0 + EA( 1 + (∂w/∂x)2 − 1) (∂w/∂x)
µ 2 +α −κ 2
= p + f (x, t).
∂t ∂t ∂t∂x ∂x 1 + (∂w/∂x)2
(2.5)
Si on considère

2 le mouvement de la corde suffisant faible par rapport à sa longueur
alors ∂w
∂x
≈ 0. Dans ce cas, on retrouve l’équation des cordes linéarisée. Soit,

∂ 2 w(x, t) ∂w(x, t) ∂ 3 w(x, t) ∂ 2 w(x, t)

µ + α − κ − T0 = f (x, t). (2.6)
∂t2 ∂t ∂t∂x2 ∂x2
Dans la suite de ce document, la discrétisation par éléments finis dans le forma-
lisme des Systèmes Hamiltoniens à Port est d’abord appliquée au cas de la corde
linéaire et est ensuite appliquée au cas non linéaire avec le modèle de corde de
Kirchhoff-Carrier.

12
3 Méthode des éléments finis
La relation (2.6) décrit une équation aux dérivées partielles (EDP). Le recours à
une méthode numérique permet d’en donner une solution approchée. La méthode
des éléments finis [13, 28], [29, §5.3], est l’une des techniques permettant cette ap-
proximation et se décompose en deux parties : formulation faible de l’EDP puis
approximation du déplacement de chaque élément du domaine par des fonctions
d’interpolation.

3.1 Formulation faible

Un problème physique est généralement décrit par un système d’équations dif-
férentielles ou d’équations aux dérivées partielles. Cette formulation est appelée
formulation forte [30]. Il est cependant possible d’exprimer ces systèmes par une for-
mulation dite faible (ou variationnelle). L’avantage de cette formulation est qu’elle
permet de se placer dans des espaces mathématiques (Sobolev, Hilbert) où des théo-
rèmes permettent de garantir l’existence et l’unicité de solutions. Se basant sur le
principe de moindre action, la formulation faible est obtenue en multipliant l’EDP
par une fonction test et en intégrant sur le domaine.
Dans la suite de cette partie, par soucis de clarté et pour une meilleure com-
préhension des sections suivantes, la méthode des éléments finis est présentée en
l’appliquant à l’équation des cordes linéarisée (2.6). Le domaine Ω = [0; L] est unidi-
mensionnel et en intégrant et en multipliant par une fonction test v, la formulation
faible de l’équation des cordes (2.6) s’exprime :
Z L 2 Z L
∂ w ∂w ∂ 3w ∂ 2w 
µ 2 +α −κ − T0 2 vdx = f vdx, (3.1)
0 ∂t ∂t ∂t∂x2 ∂x 0
L’intégration par parties des facteurs dépendant de la dérivée spatiale d’ordre deux
permet d’écrire :
Z L 2  L Z L
∂ w ∂w ∂w ∂v
2
vdx = v − dx,
0 ∂x ∂x 0 0 ∂x ∂x
Z L 3  2 L Z L 2
∂ w ∂ w ∂ w ∂v
2
vdx = v − dx.
0 ∂t∂x ∂t∂x 0 0 ∂t∂x ∂x

La fonction test v a les mêmes conditions aux limites que le problème. Dans le cas
de conditions de Dirichlet, la formulation faible de (2.6) a pour expression :

L L L L L
∂ 2w ∂ 2 w ∂v
Z Z Z Z Z
∂w ∂w ∂v
µ 2 vdx+ α vdx+ κ dx+ T0 dx = f vdx. (3.2)
0 ∂t 0 ∂t 0 ∂t∂x ∂x 0 ∂x ∂x 0

3.2 Discrétisation spatiale et interpolation polynômiale

Le principe de la méthode des éléments finis repose sur une discrétisation spatiale
du domaine en un nombre finis d’éléments. Dans le cas de la corde, le domaine Ω
est segmenté en n éléments de longueur h = Ln (figure 2).

13
 Figure 2 – Domaine divisé en n éléments de même taille

Le déplacement de chaque élément peut être approximé par interpolation polynô-

miale des déplacements des nœuds (figure 3) :

w(x, t) ≈ w1 (t)φ1 (x) + w2 (t)φ2 (x), (3.3)

où les fonctions de bases φ(x) peuvent prendre différentes expressions en fonction du
degré du polynôme choisi. En approximant la solution par des polynômes de degré
1, ces fonctions ont pour expressions :
h−x x
φ1 (x) = , φ2 (x) = . (3.4)
h h

Figure 3 – Approximation du déplacement de chaque élément

En introduisant les notations vectorielles W = [w1 w2 ]T et Φ = [φ1 φ2 ]T , le dépla-

cement en tout point s’exprime :

w(x, t) ≈ ΦT W. (3.5)
La fonction test v(x) étant arbitraire, elle peut être approximée par le vecteur conte-
nant les fonctions de bases. Soit,

v(x, t) ≈ ΦT V. (3.6)
La force extérieure f (x, t) est une fonction continue en temps et en espace et peut
aussi être approximée telle que :

f = ΦT F, (3.7)
T
où F = [F1 F2 ] , avec F1 et F2 les forces appliquées aux nœuds délimitant l’élé-
ment. La formulation faible (3.2) sur un élément s’exprime à présent :

Z h Z h Z h
T T
T T
µΦ ẄΦ Vdx + T T
αΦ ẆΦ Vdx + κΦ0 ẆΦ0 Vdx
0 0 0
Z h Z h
0T 0T
+ T0 Φ WΦ Vdx = ΦT FΦT Vdx, (3.8)
0 0

14
 Z h Z h Z h
T
T
µV ΦΦ Ẅdx + T T
αV ΦΦ Ẇdx + κVT Φ0 Φ0 Ẇdx
T
0 0 0
Z h Z h
T 0 0T
+ T0 V Φ Φ Wdx = VT ΦΦT Fdx, (3.9)
0 0

En factorisant par VT ,

Z h Z h Z h Z h
T T 0 0T 0 0T
µΦΦ dx Ẅ + αΦΦ + κΦ Φ dx Ẇ + T0 Φ Φ dx W = ΦΦT Fdx,
|0 {z } |0 {z } |0 {z } | 0 {z }
Me Ce Ke Fe
(3.10)
où Me , Ke , Ce et Fe correspondent respectivement aux matrices élémentaires de
masse, de raideur, d’amortissement et au vecteur de force. Ces variables ont pour
expressions
" Rh 2 Rh # " #
φ (x)dx φ 1 (x)φ 2 (x)dx µh 2 1
Me = µ R h 0 1 0
Rh 2 = ,
0
φ 1 (x)φ 2 (x)dx 0
φ 2 (x)dx 6 1 2
" Rh h 0
# " #
02 0
R
φ (x)dx φ (x)φ (x)dx T 0 1 −1
Ke = T0 R h 00 1 0 0
R h1 02 2 = ,
0
φ1 (x)φ 2 (x)dx 0
φ 2 (x)dx h −1 1
" #
α κ 6 F 1
Ce = Me + Ke , Fe = Me .
µ T0 µh F2

3.3 Assemblage des matrices

Pour construire les matrices, dites globales, de masse et de raideur sur l’ensemble
du domaine Ω, les matrices élémentaires doivent être assemblées. Nous présentons
ici le cas où le domaine est divisé en deux éléments (donc trois noeuds). Les matrices
élémentaires bleues correspondent à celles du premier élément et celles en rouge au
second. Le partage du second noeud par le premier et le deuxième élément conduit
à une somme des facteurs des matrices à ce point.
Les matrices globales ont alors pour expressions :
   
2 1 0 1 −1 0
µh  T0 
M=  1 2+2 1 , K =  −1 1 + 1 −1 ,
 
6 h
0 1 2 0 −1 1

     
2 1 0 1 −1 0 F1
αh   κ 6 
C=  1 2 + 2 1  +  −1 1 + 1 −1 , F = M F2 + F2 .
 
6 h µh
0 1 2 0 −1 1 F3

La formulation faible du problème peut finalement être écrite sous la forme :

MẄ + CẆ + KW = F. (3.11)

15
3.4 Conditions aux limites
Le système (3.11) représente une corde discrétisée en un nombre finis d’éléments.
Cependant, il ne prend pas en compte les conditions aux limites. La corde étant
fixée à ses extrémités, il est essentiel d’imposer des déplacements nuls aux premier
et dernier nœuds du système. Deux techniques sont alors envisagées : le recours à
des multiplicateurs de Lagrange [30, § 12.6.4] ou une modification des matrices du
système à résoudre.

3.4.1 Multiplicateurs de Lagrange

Les multiplicateurs de Lagrange, notés λi , permettent d’imposer une valeur de

déplacement à certaines nœuds et représentent alors les réactions aux appuis. Dans
notre cas, les deux conditions aux limites sont des conditions de Dirichlet. En intro-
duisant ces variables, la formulation du problème statique s’exprime :
    
K DT W F
= , (3.12)
D 0 λ ud
avec λ = [λ1 , λ2 ]T , ud = [ud1 , ud2 ]T les valeurs des déplacements imposés et D
une matrice de dimension 2 × (n + 1). Puisque la corde est fixée à ses extrémités,
ud1 = ud2 = 0 et la matrice D est remplie de valeurs nulles exceptée pour les nœuds
où le déplacement est imposée, auquel cas D(1, 1) = D(n+1, n+1) = 1. Par exemple,
si la corde est discrétisée en trois éléments (donc quatre nœuds) et en posant T0 et
h égaux à 1 dans la matrice K pour plus de clarté, le problème statique a pour
expression :

1 −1 0 0 1 0 w1
    
F1
−1 2 −1 0 0 0  w2   F2 
 0 −1 2 −1 0 0  w3   F3 
    
  =  . (3.13)
0 0 −1 1 0 1  w4   F4 

1 0 0 0 0 0   λ1  ud1 = 0
0 0 0 1 0 0 λ2 ud2 = 0
Dans le cas dynamique, la matrice M n’est pas affectée par les multiplicateurs de
Lagrange et le problème dépourvu d’amortissement s’écrit :
       
M 0 Ẅ K DT W F
+ = . (3.14)
0 0 λ D 0 λ ud
| {z } | {z } | {z } | {z }
M Ẅ K W

En supposant la solution harmonique, soit Ẅ = −ω 2 W avec ω les pulsations

propres du système, et la force extérieure nulle, la corde discrétisée peut s’écrire
comme un problème aux valeurs propres généralisé :

(−ω 2 M + K)W = 0, (3.15)

dont les valeurs propres et vecteurs propres sont respectivement les pulsations propres
au carré ω 2 et les déformées modales Λ. La figure 4 représente les déformées modales
des deux premiers modes ainsi que les réactions aux appuis afin que le déplacement
du premier et du dernier nœud soit nul.

16
 Déformées modales
0.5

Amplitude
0

Premier mode
Second mode

-0.5
0 10 20 30 40 50 60
Points discrets de la corde

Figure 4 – Deformées modales des deux premiers modes de la corde et réactions aux appuis

3.4.2 Modification des éléments du système

Un autre moyen pour introduire des conditions de Dirichlet dans le système est de
modifier la matrice de raideur K et le vecteur de force F sans augmenter le nombre
de degrés de liberté, contrairement aux multiplicateurs de Lagrange. Soit, pour un
système de taille 4 × 4 et avec T0 et h de valeur unitaire, le système
    
K11 K12 K13 K14 w1 F1
K21 K22 K23 K24  w2  F2 
K31 K32 K33 K43  w3  = F3  (3.16)
    
K41 K42 K43 K44 w4 F4
devient
    
1 0 0 0 w1 w̃1
 w2  = F2 − K21 w̃1 − K24 w̃4 
0 K22 K23 0    
 (3.17)
0 K32 K33 0 w3
   F3 − K31 w̃1 − K34 w̃4 
0 0 0 1 w4 w̃4
avec w̃1 et w̃3 les valeurs souhaitées des déplacements imposés. Les déplacements
étant nuls aux extrémités de la corde,
    
1 0 0 0 w1 0
0 K22 K23 0 w2  F2 
0 K32 K33 0 w3  = F3  (3.18)
    
0 0 0 1 w4 0
Dans le cas dynamique, il est possible d’imposer des déplacements nuls en découplant
les nœuds considérés. Pour ce faire, la matrice M est affectée de la même manière
que la matrice K. Soit,
       
1 0 0 0 ẅ1 1 0 0 0 w1 0
0 M22 M23 0 ẅ2  0 K22 K23
    0 w2  F2 
   

0 M32 M33 + = . (3.19)
0 ẅ3  0 K32 K33 0 w3  F3 
0 0 0 1 ẅ4 0 0 0 1 w4 0
La dynamique des nœuds aux conditions limites est alors
ẅi + wi = Fi , (3.20)
avec i l’indice du point considéré. Elle dépend donc des conditions initiales et des
valeurs de forces attribuées. Dans notre cas, les conditions initiales sont toujours

17
nulles quelque soit le nœud. En revanche, pour que le déplacement reste nul au
cours du temps en ces points spécifiques, il convient que fi = 0.

3.4.3 Discussion sur les deux méthodes

Les deux méthodes précédentes, permettant l’introduction des conditions aux li-
mites dans la formulation du problème, ont chacune leurs avantages et désavantages.
Les multiplicateurs de Lagrange permettent d’obtenir les valeurs des réactions aux
appuis mais le nombre de degrés de liberté du problème augmente et les propriétés
des matrices sont modifiées (la matrice M n’est plus inversible). L’autre technique a
l’avantage de ne pas augmenter le nombre de degrés de liberté mais nous n’avons pas
accès aux valeurs des réactions aux appuis et ne fonctionne que lorsque les conditions
aux limites sont nulles aux nœuds considérés.
La seconde méthode est tout de même employée dans ce document puisqu’à partir
de la partie 5, il est nécessaire d’inverser la matrice M pour modéliser la corde dans
le formalisme des SHP.

18
4 Systèmes Hamiltoniens à Ports
Cette section a pour but d’introduire le formalisme des Systèmes Hamiltoniens
à Ports [17]. L’objectif de ce formalisme est de garantir a priori que la passivité du
système est garantie au cours de la modélisation. Dans la première partie de cette
section, un cas général de formulation d’un SHP à l’aide d’une représentation d’état
est présenté. Les matrices d’interconnexion des différents composants (stockants,
dissipatifs et de ports) sont identifiées et la conservation de l’énergie est assurée
à partir de ces éléments. La seconde partie présente les notions géométriques sur
lesquels le formalisme des SHP s’appuie. Les notions de flux et d’efforts sont ainsi
abordées au travers de la structure de Dirac. La troisième partie applique les aspects
théoriques à un exemple simple qu’est l’oscillateur mécanique amorti.

4.1 Cas général

Le formalisme des SHP est basé sur l’interconnection d’éléments stockants, dis-
sipatifs et de ports. Cette connexion est assurée par une matrice S, avec S = −S T .
La formulation dite algébro-différentielle des SHP [31] s’écrit :
 dz   
dt
∇H
 d  = S  z(d)  (4.1)
−y u
La variable d’état z (éléments stockants) est ainsi liée au gradient de l’énergie du
système H (∇ correspond à l’opérateur gradient), l’amortissement d à une fonction
z(d) et la sortie y à une entrée u (le produit y T u quantifiant la puissance apportée
au système).
Une autre formulation des SHP, gouvernée par la représentation d’état [17] peut être
obtenue à partir du système (4.1). Soit,
(
dz


dt
= J(z) − R(z) ∇H(z) + G(z)u
(4.2)
y = G(z)T ∇H(z)
où la matrice J = −J T est anti-symétrique, R = RT ≥ 0 est une matrice positive
et G permet de choisir sur quelle variable d’état la force est appliquée. La variation
d’énergie de ce système au cours du temps s’exprime donc :

dE dz
= ∇H(z)T = ∇H(z)T J(z)∇H(z) − ∇H(z)T R(z)∇H(z) + y T u , (4.3)
dt dt | {z } | {z } |{z}
Pc =0 Q≥0 Pext

Pc traduit les échanges de puissance entre les composants stockants et est nulle
puisque J est anti-symétrique. Q représente la puissance dissipée et est positive
ou nulle puisque R ≥ 0. Le dernier terme Pext correspond à la contribution de
la puissance extérieure au système. Finalement, la variation d’énergie du système
est soit nulle (cas conservatif), soit négative (cas dissipatif), mais jamais positive,
puisqu’aucun système ne crée de l’énergie spontanément.

4.2 Structure de Dirac

Le formalisme des Systèmes Hamiltoniens à Ports s’appuie sur la notion de struc-
ture de Dirac ([17, §2.2],[32, §2.1.2], [33]), dont l’une des propriétés fondamentale est

19
la conservation de puissance. Cette structure permet de mettre en relation des flux
f et des efforts e de sorte que la puissance totale eT f soit nulle. En notant F l’espace
vectoriel contenant les éléments caractérisant les flux f ∈ F et E := F ∗ l’espace dual
de F contenant les efforts e ∈ E, il est possible de définir l’espace des variables de
ports F × E. Il existe une forme bilinéaire << >> définie par :

<< (fa , ea ), (fb , eb ) >>:= hea |fb i + heb |fa i, (4.4)

avec (fa , ea ), (fb , eb ) ∈ F × E et h|i le produit dual. Finalement, une structure de
Dirac sur F × E est un sous-espace D ⊂ F × E tel que :

D = D⊥ , (4.5)
où D⊥ représente l’ensemble des éléments de F × E orthogonaux à l’ensemble D.
La variation d’énergie (4.3) du système (4.2) peut à présent être réécrite sous la
forme :

Ḣ(z) = ∇H(z)T ż(t) = h∇H(z)|ż(t)i, (4.6)

où ż est défini comme le vecteur des variables de flux f et ∇H(z) comme le vecteur
des variables d’efforts e. Avec ces nouvelles notations, la représentation d’état (4.2)
s’écrit :
(
f = (J − R)e + Gu,
(4.7)
y = GT e

4.3 Exemple : Oscillateur mécanique amorti

Le système considéré, représenté en figure 5, est un oscillateur mécanique amorti
constitué d’un ressort de raideur k, d’une masse m, d’un amortisseur c, excité par
une force extérieure u.

Figure 5 – Système masse-ressort-amortisseur linéaire soumis à une force extérieure

En reprenant les notations établies dans la représentation d’état (4.2), z = [q, p]T , où
q et p = m dq
dt
correspondent respectivement à l’élongation du ressort et à la quantité
de mouvement de la masse. L’énergie H(z) a pour expression H = 21 z T Bz, avec
B = diag(k, m−1 ). L’énergie potentielle du ressort s’exprime, selon la loi de Hooke,
en fonction de l’élongation q du ressort. Dans le cas linéaire,
(
q̇ = −fk ,
d 1 2
 (4.8)
ek = dq 2 kq ,
avec −fk la vitesse du ressort au point d’attache avec la masse (le signe − vient du
fait qu’on veut que le produit fk ek soit une puissance entrante) et ek la force de

20
rappel du ressort à ce même point. L’énergie cinétique s’exprime en fonction de la
quantité de mouvement de la masse. Soit,
(
ṗ = −fm ,
d 1 2
 (4.9)
em = dp 2m
p ,
avec −fm la force exercée sur la masse et où em = mp correspond à la vitesse de la
masse. Les systèmes masse et ressort peuvent être couplés en appliquant le principe
d’action-réaction de Newton :
(
−fk = em ,
(4.10)
fm = ek .
Le système sans amortissement ni force extérieure s’exprime :
     ∂H 
q̇ 0 1 ∂q
= (4.11)
ṗ −1 0 ∂H ∂p
|{z} | {z } | {z }
f J e

Selon [31], en présence d’amortissement et excité par une force extérieure, les ma-
trices du système (4.2) ∇H(z), J(z), R(z) et G(z) ont pour expressions :

       
kq 0 1 0 0 0
∇H(z) = , J= , R= , G= .
m−1 p −1 0 0 c 1
La simulation de ce système par une méthode numérique adaptée (cf. section 5.4,
[31]) est représentée en figure 6 avec pour paramètres m = 0.1 kg, k = 5 N, c = 0.1
N.m/s. Les conditions initiales sont z = [0, l0 ] avec l0 = 0.1 m. Il y a absence de
force extérieure hormis sur l’intervalle [2s, 5s] où Fext = kl0 /2 et sur l’intervalle [10s,
12s] où Fext = kl0 .

0.02

0.01

-0.01
0 5 10 15
t(s)
0.1

-0.1
0 5 10 15
t(s)

100

10-5

10-10
0 5 10 15
t(s)

Figure 6 – Simulation d’un masse-ressort linéaire

21
5 Modélisation d’une corde linéaire par la méthode des élé-
ment finis dans le formalisme des Systèmes Hamiltoniens
à Ports
Après avoir détaillé dans les parties précédentes la méthode des éléments finis et
le formalisme des Systèmes Hamiltoniens à Ports séparément, cette partie montre
comment il est possible de discrétiser spatialement un système physique (une corde
dans le cas présent) par la méthode des éléments finis et de l’exprimer dans le
formalisme des SHP. La première partie de cette section montre comment l’énergie
totale (l’hamiltonien) du système peut être formulée à partir de l’approximation par
éléments finis. Le modèle de corde est ensuite écrit dans le formalisme des Systèmes
à Ports. La dernière partie détaille la méthode numérique garantissant la passivité
du système en temps discret.

5.1 Formulation faible du modèle linéaire

La formulation faible du problème linéaire a été réalisée dans la partie 3 (équation
(3.11)) et nous en rappelons ici le résultat. En multipliant l’équation des cordes
linéarisée par une fonction test v et en intégrant sur le domaine Ω, le formulation
faible a pour expression :

MẄ + CẆ + KW = F. (5.1)

5.2 Approximation des énergies

La dynamique d’une corde vibrante non amortie et sans force extérieure est régie
par l’équation d’onde

∂ 2w ∂ 2w
µ − T 0 = 0. (5.2)
∂x2 ∂t2
L’énergie totale de la corde est la somme des énergies cinétique et potentielle.
Soit ([34, §4.2.1]),
!
1 L
Z  ∂w 2  ∂w 2
H= µ + T0 dx. (5.3)
2 0 ∂t ∂x
En approximant le déplacement par w(x, t) = ΦW, l’énergie totale s’exprime

1 L
Z
H= µ(ΦẆ)2 + T0 (WΦ0 )2 dx, (5.4)
2 0
1 L
Z
T
H= µẆT ΦT ΦẆ + T0 WT Φ0 Φ0 Wdx,
2 0
Z L Z L !
1 T
H= ẆT µΦΦT dx Ẇ + WT T0 Φ0 Φ0 dxW ,
2 0 0

L’énergie totale du système s’exprime finalement :

1 T T

H= Ẇ MẆ + W KW (5.5)
2

22
5.3 Ecriture dans le formalisme des SHP
Le but est à présent de pouvoir lier la formulation faible du problème à son énergie
H et de se rapprocher du formalisme hamiltonien en termes de flux et d’efforts. En
notant l’état du système z = [W Ẇ], l’entrée u = F, la sortie y = Ẇ et en
définissant le vecteur des variables de flux f = dz
dt
et le vecteur des variables d’efforts
T
e = ∇H = [KW MẆ] , le système peut s’écrire dans la même configuration que
(4.2). Soit,

" # " # " #!" # " #

−1
 Ẇ 0 M 0 0 KW 0
= − +

 F
−M−1 0 M−1 CM−1 M−1 |{z}


 Ẅ 0 MẆ
} | {z } | {z } u



 | {z } | {z } | {z
 f=ż "J R e=∇H G

#
h i KW (5.6)
 Ẇ = 0 M−1
| {z } MẆ

 |{z}


 y
GT

 | {z }


 e=∇H
+ Conditions aux limites comprises dans M et K (cf. 3.4)

Le système (5.6) représente ainsi la corde vibrante discrétisée par éléments finis
comme un Système Hamiltonien à Ports.

5.4 Méthode numérique à passivité garantie

L’intérêt principal des SHP est de garantir la passivité du système, c’est à dire
qu’il n’y a pas création d’énergie spontanée. Pour simuler numériquement le système
dynamique (5.6), une version discrète du gradient de l’hamiltonien doit être formulée.
Une approche possible est de s’appuyer sur un schéma d’Euler ([31], [22]) restituant
une version discrète de la règle de la dérivation en chaîne pour la fonction composée
E = H ◦ z. Soit,
δE(t, δt)
T δz(t, δt)
= ∇d H z, δz(t, δt) , (5.7)
δt δt


avec δt le pas temporel, ∇d H z, δz(t, δt) le gradient discret du hamiltonien H et
où δz(t,δt)
δt
approxime le vecteur z(t) par différenciation décentrée avancée, tel que
δz(t, δt) z(t + δt) − z(t)
= . (5.8)
δt δt
Dans le cas d’un hamiltonien à variables séparées de la forme H(z) = N
P
n=1 Hn (zn ),
le gradient discret associé s’exprime
(
Hn (zn +δzn )−Hn (zn )
δzn
si δzn 6= 0,
[∇d H(z, δz)]n = 0
(5.9)
Hn (zn ) sinon.
    
K 0 W K 0

En utilisant cette définition avec H = = z , le gradient
0 M Ẇ 0 M
discret d’un hamiltonien quadratique régénère la méthode du point milieu :
1 1
K 0 2 (z + δz)2 − 2 z 2

∇d H(z, δz) = .
0 M δz

23
  
K 0
En posant B = ,
0 M
1 2
z + zδz + 21 δz 2 − 12 z 2
∇d H(z, δz) = B 2 ,
δz
 
 δz  q + δq/2
∇d H(z, δz) = B z + =B . (5.10)
2 p + δp/2
En remplaçant ∇d H(z, δz) par son expression (5.10) dans la première équation du
système (4.2), on obtient :
δz  δz 
= (J − R)B z + + Gu,
δt 2
h  δz  i
δz = (J − R)B z + + Gu δt,
2
δt  
δz − (J − R)Bδz = (J − R)Bz + Gu δt.
2
A chaque instant k, l’expression de la variation δz est de la forme :
δt −1  
δz(k) = I − D Dz(k) + Gu(k) δt, (5.11)
2
avec D = (J − R)B. La résolution est explicitée dans l’algorithme 1.

Algorithme 1 :
Données : J, R, B, G, x0 , uk
Résultats : z(k), y(k)
1 pour k = 0 à kt faire
2 calculer δz(k) avec l’équation (5.11);
3 z(k + 1) ← z(k) + δz(k) ;
4 y(k) ← GT ∇Hd (z(k)) ;

5.5 Simulation numérique

Une fois la résolution numérique effectuée en employant le gradient discret décrit
en section 5.4, l’état z(k) et la sortie y(k) de chaque nœud sont déterminés pour
chaque instant k. A titre d’exemple, la représentation temps-fréquence de la sortie
du 3ème nœud d’une corde excitée par une force ponctuelle extérieure (détaillée en
Annexe 9.1) est affichée en figure 7a. La corde,
q discrétisée en 60 éléments, d’une
1 T0
longueur de L = 0.6 m, est accordée à f0 = 2L µ
= 110 Hz et pincée au cinquième
de sa longueur. Les valeurs des amortissements fluide et structurel sont choisies
empiriquement [35] tels que α = 3µ kg.m-1.s-1 et κ = 0.01µ kg.m.s-1. La figure 7b
représente les énergies du système, et on remarque que l’énergie totale décroît au
cours du temps dû à l’amortissement de la corde.
Les figures 8a et 8b représentent le cas où la corde est dépourvue d’amortissement,
soit α = κ = 0. L’énergie totale H, représentée en figure 8b, est constante dans ce cas,
prouvant ainsi la stabilité de la simulation qui ne crée pas d’énergie spontanément.

24
 10-4
4

3.5

3
Energie potentielle
2.5 Energie cinetique

Energie (J)
Energie totale
2

1.5

0.5

0
0 0.5 1 1.5 2 2.5 3
Temps (s)

(a) Représentation temps-fréquence de la vitesse du

(b) Energies du système
noeud n◦ 3

Figure 7 – Cas avec un amortissement C non nul

10-4
4

3.5

2.5
Energie (J)

2 Energie potentielle
Energie cinetique
1.5
Energie totale
1

0.5

0
0 0.5 1 1.5 2 2.5 3
Temps (s)

(a) Représentation temps-fréquence de la vitesse du

(b) Energies du système
noeud n◦ 3

Figure 8 – Cas avec un amortissement C nul

25
6 Extension au cas de la corde non linéaire
6.1 Modèle de Kirchhoff-Carrier
Le modèle linéaire de corde vibrante étudié dans les sections précédentes propose
une bonne approximation du comportement du système tant que les amplitudes de
déplacement sont faibles. Lorsque celles-ci deviennent grandes, ce modèle n’est plus
suffisant, ne prenant pas en compte les variations de tension dues aux variations
de longueur de la corde. A partir de l’équation des cordes (2.4) développée dans la
section 2, il est possible de prendre en compte ce phénomène. Une autre écriture
de l’équation (2.4) est possible en prenant en compte la non-linéarité de manière
globale (et non locale comme c’est le cas dans l’équation (2.4)) et s’exprime [36] :
Tension globale
z }| !{
2 Z L
∂ w EA ∂w 2 ∂ 2w
µ = T0 + dx , (6.1)
∂t2 2L 0 ∂x ∂x2
| {z }
terme non linéaire

La tension totale est constituée de la tension statique T0 (équivalente au cas linéaire)

sommée à un terme non linéaire qui traduit la variation de tension dans la corde.

6.2 Formulation faible du modèle non linéaire

La formulation faible du problème linéaire s’obtient de la même manière que pour
le cas linéaire, à savoir en multipliant par une fonction v test et en intégrant sur le
domaine :
Z L" 2 Z L  2 ! #
∂ w ∂w ∂ 2w
µ 2 − 1+ dx T0 2 vdx = 0, (6.2)
0 ∂t 0 ∂x ∂x

EA
R L  ∂w 2
avec  = 2LT 0
. En notant T (t) = 1+ 0 ∂x
dx, l’intégration par parties du terme
dépendant de la dérivée spatiale d’ordre deux a pour expression :
Z L Z L
∂ 2w h ∂w iL ∂w ∂v
T (t) T0 2 vdx = T (t) T0 v − T (t) T0 dx. (6.3)
0 ∂x ∂x 0 0 ∂x ∂x
v satisfaisant les conditions aux limites, soit v(0) = v(L) = 0, la formulation faible
s’exprime finalement :
Z L 2 Z L
∂ w ∂w ∂v
µ 2 vdx + T (t) T0 dx = 0. (6.4)
0 ∂t 0 ∂x ∂x
Dans la section 3.2, il a été montré que le déplacement w(x, t) pouvait être
approximé par interpolation polynomiale à l’aide de fonctions de bases tel que
w(x, t) ≈ ΦT W et que v ≈ ΦT V. En utilisant ces approximations, les termes
linéaires peuvent se mettre sous les formes :
Z L 2 Z L
∂ w T
µ 2 vdx = V µΦΦT dxẄ = VT MẄ, (6.5)
0 ∂t 0
Z L Z L
∂w ∂v T
T0 dx = V T
T0 Φ0 Φ0 dxW = VT KW. (6.6)
0 ∂x ∂x 0

26
Le terme non linéaire T (t), sous les mêmes approximations, s’exprime [37] :
Z L  2 Z L
∂w T
T (t) = 1 +  dx = 1 +  (Φ0 W)2 dx, (6.7)
0 ∂x 0
Z L
T
T (t) = 1 +  WT Φ0 Φ0 Wdx,
0
Z L
 T
T (t) = 1 + WT T0 Φ0 Φ0 dxW,
T0 0

T (t) = 1 + WT KW, (6.8)
T0
où les matrices M et K sont les matrices globales issues de l’assemblage des ma-
trices élémentaires Me et Ke définies à la section 3.2. La méthode des éléments finis
appliquée au modèle de Kirchhoff-Carrier permet d’écrire :
  
VT MẄ + VT 1 + WT KW KW = 0. (6.9)
T0
En factorisant par VT , en introduisant le facteur d’amortissement et la force exté-
rieure de la même manière que dans le cas linéaire et en posant β = T0 , le modèle
non linéaire s’exprime finalement :
 
T
MẄ + CẆ + 1 + βW KW KW = F. (6.10)

6.3 Approximation des énergies du modèle non linéaire

Les énergies cinétique et potentielle d’une corde vibrante modélisée par le modèle
de Kirchhoff-Carrier ont pour expressions :

1 L  ∂w 2
Z
Ec = µ dx, (6.11)
2 0 ∂t
Z L  2 ! Z L  2
∂w 1 ∂w
EpN L = 1 +  dx T0 dx. (6.12)
0 ∂x 2 0 ∂x
L’énergie cinétique est identique au cas linéaire. Dans le cas non linéaire, seule l’éner-
gie potentielle est modifiée en raison du facteur de tension variable. En approximant
le déplacement par w(x, t) ≈ ΦT W, l’énergie cinétique a pour expression
1
Ec = ẆT MẆ, (6.13)
2
et l’énergie potentielle non quadratique s’exprime
Z L ! Z
T 1 L T
EpN L = 1 +  (Φ0 W)2 dx T0 (Φ0 W)2 dx, (6.14)
0 2 0
Z L !
1 L
Z
T 0 0T T
EpN L = 1 +  W Φ Φ Wdx T0 WT Φ0 Φ0 Wdx,
0 2 0
!
 1
EpN L = 1 + WT KW WT KW. (6.15)
T0 2

27
L’énergie totale du modèle non linéaire étant la somme des énergies cinétique et
potentielle :
!
1 T

T
1
T
HN L = Ẇ MẆ + 1 + βW KW W KW . (6.16)
2 2

6.4 Gradient discret de l’hamiltonien non quadratique

Contrairement au cas linéaire, l’hamiltonien du problème non linéaire défini à
l’équation (6.16) n’est pas à variables séparées. L’expression de la version discrète
du gradient de l’hamiltonien proposée à l’équation (5.9) n’est donc plus valable. D’un
point de vue général, la version discrète de ∇H(z)T ż doit être égale à la variation
d’énergie entre deux instants successifs ([38], [39]). Soit,

∇d H(z, δz)T δz = H(z + δz) − H(z), (6.17)

avec ∇d H(z, δz) la version discrète du gradient de l’hamiltonien. Pour un système
de dimension N , la variation d’énergie peut se décomposer comme N sommes tel
que

H(z + δz) − H(z) = H(z1 + δz1 , z2 + δz2 , ..., zN + δzN ) − H(z1 , z2 + δz2 , ...)

+H(z1 , z2 + δz2 , ...) − H(z1 , z2 , z3 + δz3 , ...)

.
+...

+H(z1 , ..., zN −1 , zN + δzN ) − H(z1 , z2 , ..., zN )

(6.18)
Cette formulation dépend cependant de l’ordre d’indexation χ des variables d’état
de tel sorte que
 
H(z1 + δz1 ) − H(z1 , z2 + δz2 , ...)
δx1
[∇d H(z, δz)]χ =  ... . (6.19)
 
H(z1 , ..., zN −1 , zN + δzN ) − H(z1 , ..., zN −1 , zN )
δzN

Finalement, une expression unique du gradient discret est retenue en moyennant la

somme des chemins par le nombre de permutations possibles. Soit,
1 X
∇d H(z, δz) = [∇d H(z, δz)]χ (6.20)
N! χ
Cette définition du gradient discret est ainsi valable dans le cas linéaire et non
linéaire. (Il a été prouvé dans [39, A1], que ce gradient discret régénérait la méthode
du point milieu lorsque le Hamiltonien est quadratique).

28
6.5 Résolution et simulation numérique
La représentation d’état discrète du système à résoudre pour chaque instant k
est de la forme :
(
dz(k)


dt
= J − R ∇Hd z(k), δz(k) + Gu(k)

(6.21)
y(k) = G(z)T ∇H z(k), δz(k)
Par rapport au cas linéaire, seul l’hamiltonien est modifié par une tension addition-
nelle variable dans le temps comme expliqué à la section précédente. Pour rappel,
les autres matrices ont pour expressions :

M−1
     
0 0 0 0
J= , R= , G= , u = F, (6.22)
−M−1 0 0 M−1 CM−1 M−1
où J est bien antisymétrique et R est positive ou nulle.
La définition du gradient discret émise à l’équation (6.20) n’est pas aisée à cal-
culer, surtout dans un but d’une synthèse sonore en temps réel, puisqu’elle dépend
directement du nombre d’éléments employés pour discrétiser la corde qui peut être
élevé. Afin de pallier à ce problème, il est possible de trouver une nouvelle expres-
sion de ce gradient en profitant de la structure de l’énergie potentielle. En notant
V = WT KW l’énergie potentielle du cas linéaire, celle du cas non linéaire peut
s’écrire

EpN L = F ◦ V (W), avec F (V ) = V + βV 2 . (6.23)

En appliquant la règle de dérivation en chaîne au gradient de l’hamiltonien, celui-
ci s’exprime [22] :
 
∇d F (V, δV )∇d V (W, δW)
∇d H(z, δz) = , (6.24)
MẆ
avec δV = ∇d V (W, δW)T δW, ∇d F (V, δV ) = 1 + 2β(V + δV /2) et ∇d V (W, δW) =
K(W + δW/2)
La méthode de résolution consiste à trouver le point fixe δz en appliquant une
méthode itérative (d’indice i) à chaque instant k permettant de trouver la modula-
tion de tension comprise dans le facteur ∇d F (V, δV ). L’algorithme pour chaque pas
de temps est le suivant [22] :
— A. Calculer F(i = 0) = ∇d F (V, 0) ;
— B(i). Calculer δz (i) avec (6.21) dans lequel ∇d F (V, δV ) est remplacé par F(i) ;
— C(i). Calculer δV (i) = ∇d V (W, δW(i) )T δW(i) ;
— D(i). Calculer F(i + 1) = ∇d F (V, δV (i) ) ;
— E(i). Si F(i + 1) − F(i) < seuil, la première équation du système (6.21) est
approximée par δz = δz (i) , sinon on incrémente i et on reprend à l’étape B(i).

29
6.6 Résultats
Les simulations sont réalisées avec une corde discrétisée en 40 éléments, fixée à ses
extrémités, de longueur L = 0.65 m, de diamètre d = 0.79 mm, de masse volumique
ρ = 7800 kg.m-3, de module d’Young E = 210 GPa. La corde est excitée ponctuel-
lement au cinquième de sa longueur par une force excitatrice simulant l’action du
musicien, c’est à dire une rampe jusqu’à une valeur Fmax sur une durée de 10 ms.
La valeur de Fmax est modifiée en fonction de la hauteur de lâcher choisie (cf 9.1) .
Les représentations temps-fréquences de la vitesse d’un des nœuds affichées sur les
figures 9a, 9b et 10a permettent d’observer les effets non linéaires dus aux différentes
hauteurs de lâcher. La modulation de tension crée un pitch glide effect qui se traduit
par une augmentation des fréquences propres de la corde lorsque son déplacement
est grand. Les fréquences tendent vers les valeurs du cas linéaire lorsque l’énergie de
la corde diminue en raison de l’amortissement. De plus, un enrichissement spectral
est présent. Le cas linéaire prévoit que l’harmonique multiple du point de pince-
ment (le cinquième harmonique ici) est absent. Dans le cas non linéaire, l’énergie
de cette harmonique est d’autant plus grande que le déplacement est important.
La figure 10b représente les énergies du système. La corde étant amortie, l’énergie
totale décroît au cours du temps.

(a) (b)

Figure 9 – Spectrogramme de la vitesse de la corde pour différentes hauteurs de lâcher

0.07

0.06

0.05 Energie potentielle

Energie cinétique
Energie (J)

0.04
Energie totale
0.03

0.02

0.01

0
0 0.5 1 1.5 2 2.5 3
Temps (s)

(a) (b)

Figure 10 – A gauche, spectrogramme de la vitesse de la corde pour h = 0.01 m. A droite, les

énergies du système lorsque h = 0.01 m.

30
La figure 11 correspond aux énergies cinétique, potentielle et totale de la corde
lorsque celle-ci est dépourvue d’amortissement et lâchée à une hauteur h = 0.01 m.
L’énergie totale est constante, prouvant ainsi que la méthode de résolution garantit
aussi la passivité du système lorsque celui-ci est non linéaire.

0.07

0.06

0.05
Energie (J)

0.04 Energie potentielle

Energie cinétique
0.03
Energie totale
0.02

0.01

0
0 0.5 1 1.5 2 2.5 3
Temps (s)

Figure 11 – Energies du modèle de corde non linéaire avec un amortissement nul

31
7 Vers un démonstrateur en temps-réel
La synthèse sonore par modèle physique permet de modéliser des objets d’une
manière fine mais au prix d’un temps de calcul important. A l’heure actuelle, peu
d’instruments virtuels sont modélisés dans ce paradigme. Modalys et les produits
développés par les sociétés AAS (Applied Acoustics Systems, [40]), MODARTT
(Models and Data for Arts and Technology, [41]) ou encore ARTURIA [42] font
alors figure d’exception.
En vue d’une future intégration de résonateurs non linéaires dans Modalys, un
démonstrateur en temps-réel de la corde de Kirchhoff-Carrier est envisagé. Il s’agit
donc de pouvoir traduire le script MATLAB [43] dans un langage compilé. Modalys
étant codé en C++, notre choix se positionne sur ce langage. De plus, afin de fournir
à l’utilisateur un outil simple de démonstration, le code C++ est intégré au frame-
work JUCE pour permettre de contrôler les paramètres du modèle via une interface
graphique.

7.1 Le framework JUCE

L’environnement JUCE [44], créé par Julian Storer en 2004, est un framework
multi-plateforme écrit en C++ permettant de concevoir des applications audio
(VST, AU, AAX, RTAS, Standalone). De plus, la bibliothèque JUCE dispose de
nombreuses classes permettant la création rapide d’interfaces graphiques.

7.1.1 Application audio Standalone

Dans le cadre de ce projet, afin de fournir un démonstrateur en temps-réel, nous

créons à partir de la bibliothèque JUCE une application audio Standalone. A la
différence des plug-ins audio (VST, AU, AAX, RTAS) utilisés au sein de logiciels
host tels que les DAW (Digital Audio Workstation), une application Standalone
est un exécutable ne nécessitant aucun autre logiciel pour fonctionner. Les modèles
JUCE d’application audio Standalone sont initalement composés de trois fichiers :
— MainComponent.h
— MainComponent.cpp
— Main.cpp
L’implémentation de notre projet se réalise dans les deux fichiers MainComponent.
Le header MainComponent.h contient les attributs et les prototypes de la classe
MainComponent et le fichier source MainComponent.cpp contient la définition des mé-
thodes et leur implémentation. MainComponent hérite de la classe AudioAppComponent
qui possède par défaut trois fonctions virtuelles pures :
— AudioAppComponent::prepareToPlay()
— AudioAppComponent::releaseResources()
— AudioAppComponent::getNextAudioBlock()
Dans le cas d’une simulation en temps-différé, le temps de simulation est initialisé
dans le script et l’ensemble des calculs est effectués sur cette période. A l’inverse,
dans le cas d’une application temps-réel, il est inenvisageable de réaliser l’ensemble
des opérations avant de les envoyer à l’utilisateur. Lorsque le signal audio est généré
à une certaine fréquence d’échantillonnage (44.1 kHz par exemple), il est stocké au
fur et à mesure dans des blocs appelés buffers d’une taille spécifique (512 échantillons

32
par exemple). A chaque fois que le buffer est rempli, ses valeurs sont transmises à
la carte son qui les envoie à l’utilisateur. La fonction getNextAudioBlock() per-
met de continuer ce processus tant que l’application est ouverte. A chaque fois que
le buffer a transmis ses valeurs, le fonction regarde s’il y a encore des valeurs qui
lui parviennent, et si oui, le buffer se remplit de nouveau avec les échantillons en-
trants. La fonction prepareToPlay() est appelé au début de la simulation avant que
getNextAudioBlock() ne soit appelée et permet d’initialiser les variables déclarées
dans MainComponent.h. Pour finir, la fonction releaseResources() est appelée
lorsqu’il n’y plus d’échantillons à traiter (donc lorsque getNextAudioBlock() n’est
plus appelée) et libère la mémoire allouée lors de la simulation.
L’un des avantages de JUCE est de permettre la conception d’interfaces gra-
phiques. Dans le cas d’une application audio, les éléments de l’interface sont aussi
implémentés dans la classe MainComponent. Ils sont déclarés dans MainComponent.h
(par exemple pour un curseur : Slider nomDuSlider) et initialisés dans le construc-
teur MainComponent() situé dans le fichier MainComponent.cpp. A titre d’exemple,
un curseur réalisé avec JUCE est présenté en figure 12.

Figure 12 – Curseur réalisé avec JUCE

Les éléments graphiques de JUCE nous permettent alors de pouvoir modifier les
valeurs des paramètres que l’on veut contrôler en temps-réel sans avoir à recompiler
le projet pour un nouveau jeu de valeurs.

7.2 Eigen : une bibliothèque d’algèbre linéaire

Pour rappel, dans le cas linéaire ou non linéaire, le système à résoudre est de la
forme :
(
dz(k)



dt
= J − R ∇Hd z(k), δz(k) + Gu(k)

(7.1)
y(k) = G(z)T ∇H z(k), δz(k) ,
où l’ensemble des paramètres sont des vecteurs ou des matrices. La conversion du
script MATLAB en C++ nécessite donc de pouvoir écrire et réaliser des opérations
d’algèbre linéaire dans ce langage. Après une étude de différentes bibliothèques open-
source d’algèbre linéaire, la bibliothèque Eigen [45] s’avère être la plus adaptée tant
du point de vue de la performance (cf. Annexe 9.2) que du point de vue de la sim-
plicité d’adaptation par rapport au script MATLAB. En effet, la documentation
de cette bibliothèque fournit une transcription entre les fonctions propres à MAT-
LAB et celles propres à Eigen [46]. Le tableau 1 ci-dessous renseigne certaines des
transcriptions qui nous sont le plus utiles dans ce projet.

33
 Fonctions MATLAB Eigen
Création d’un vecteur nul A de
A = zeros(n,1) A = VectorXd::Zero(n)
dimension n
Création d’une matrice nulle A
A = zeros(n,m) A = MatrixXd::Zero(n,m)
de dimension n × m
Matrice C égale au produit
C=A∗B C=A∗B
matricielle de 2 matrices A et B
Inverse de la matrice C inv(C) C.inverse()

Tableau 1 – Transcriptions MATLAB/Eigen

7.3 Implémentations et résultats

Les implémentations des modèles de corde linéaire et non linéaire sont effectuées
sur un Macbook Pro 15 pouces équipé d’un processeur Intel Core i7 cadencé à 2.9
GHz (Turbo Boost jusqu’à 3.8 GHz) et de 16 Go de mémoire vive LPDDR3 à 2133
MHz.
Comme expliqué dans la section 7.1.1, la classe MainComponent est dôtée de trois
méthodes virtuelles pures. L’ensemble des variables sont initialisées dans la fonction
prepareToPlay() et les algorithmes de résolution temporelle sont implémentés dans
la fonction getNextAudioBlock(). La fréquence d’échantillonnage est par défaut
égale à 44.1 kHz et la taille du buffer est de 512 échantillons.
Le modèle linéaire ne nécessitant que peu de calculs à chaque pas de temps, son
utilisation en temps-réel est possible. Pour une corde discrétisée en 40 éléments, le
CPU atteint 60% d’activité. D’un autre côté, la méthode numérique employée dans
le modèle non linéaire nécessite beaucoup plus de calculs à chaque pas de temps
puisque qu’il est nécessaire de calculer par itération la modulation de tension. De
ce fait, le temps-réel n’est possible que pour un nombre très restreint d’éléments
(moins de 5) et le son de corde perçu n’est alors pas réaliste.

34
8 Conclusion et perspectives
Le logiciel de synthèse sonore Modalys permet de concevoir des instruments vir-
tuels à partir de résonateurs linéaires. Afin d’étendre les possibilités offertes par ce
logiciel et accroître le réalisme sonore des objets utilisés, l’intégration de résonateurs
non linéaires est envisagée. Une investigation de cette intégration a été réalisée au
cours de ce stage. Le problème majeur des systèmes non linéaires est leur instabilité
numérique due aux méthodes de résolution temporelle employées. Afin de pallier
à ce problème, les travaux de l’équipe S3AM s’orientent depuis quelques années
vers le formalisme des Systèmes Hamiltoniens à Ports. Reposant sur une approche
énergétique, ce formalisme permet de garantir le bilan de puissance à temps dis-
cret de systèmes non linéaires et ainsi assurer a priori la stabilité de la simulation
numérique.
Ce stage a permis de rappeler l’établissement du modèle non-linéaire de corde
vibrante. A partir de ce modèle, la méthode des éléments finis a été présentée dans
un premier temps pour le cas linéarisé conduisant à un modèle de dimension finie. De
plus, les énergies du systèmes ont elles aussi été discrétisées à partir de cette méthode.
Le formalisme des Systèmes Hamiltoniens à Ports a permis de lier l’énergie totale du
système à la formulation établie par la méthode des éléments finis. L’introduction
du gradient discret a rendu possible la mise en place d’une méthode numérique
garantissant la passivité du système linéaire. L’approximation par la méthode des
éléments finis a ensuite été appliquée de la même manière au modèle de corde de
Kirchhoff-Carrier. L’écriture du modèle non linéaire diffère par l’expression de son
énergie potentielle non quadratique traduisant une modulation de tension dans la
corde lorsque les déplacements sont grands. Grâce à la mise en place d’une seconde
méthode numérique, la simulation du modèle de corde non linéaire a été rendue
possible et sa stabilité assurée. En vue d’une application temps-réel, les modèles
ont été transcrits dans le langage C++. Cependant, même si le modèle linéaire est
utilisable en temps-réel, le modèle non-linéaire est trop coûteux en temps de calcul
dès que le nombre d’éléments augmente.
A l’issue de ce stage, certaines perspectives sont envisageables. La modélisation
d’objets non linéaires de dimensions supérieures (2D, 3D) tels que les membranes ou
les poutres pourrait être mise en place de manière identique à ce qui a été fait pour
le modèle de corde en une dimension. Dans le cas d’objets à géométries complexes, il
serait intéressant d’employer les logiciels standards d’éléments finis (ANSYS, Free-
FEM, ...) pour en extraire les matrices de masse et de raideur. Il serait aussi in-
téressant d’imposer des conditions aux limites quelconques et de pouvoir obtenir
les réactions aux limites. Cette amélioration permettrait ainsi d’interconnecter plu-
sieurs sous-structures pour en obtenir une globale. Un autre point à explorer est la
mise en œuvre d’une nouvelle méthode numérique dans le cas non linéaire car celle-
ci est trop coûteuse en temps de calcul. De plus, une optimisation des calculs est
possible. En effet, certaines des matrices comportent de nombreuses valeurs nulles
et de nombreux calculs réalisés ne sont alors pas nécessaires.

35
9 Annexes
9.1 Annexe 1 : Force extérieure
Modèle linéaire

Pour rappel, après discrétisation par la méthode des éléments finis, le système
linéaire s’exprime :

MẄ + CẆ + KW = F. (9.1)

Le vecteur de force F représente les efforts extérieurs appliqués à chacun des nœuds
du système discrétisé. Lorsque le musicien apporte de l’énergie à la corde, il tire
sur celle-ci avant de la relâcher d’une certaine hauteur. En considérant le plectre ou
le doigt comme un point (ponctuel en espace), l’effort s’applique seulement sur un
seul nœud. Pour simuler l’effort exercé par le doigt sur la corde, on considère que la
valeur de la force à ce nœud évolue linéairement jusqu’à une valeur maximale Fmax
pendant une durée de 10 ms et est nulle après le lâcher (figure 13).

0.8

0.6
Force (N)

0.4

0.2

0
0 0.005 0.01 0.015 0.02
Temps (s)

Figure 13 – Évolution de la force au point de pincement

Plutôt que d’attribuer une valeur maximale de force peu représentative pour l’uti-
lisateur, il est possible de calculer cette valeur en fonction de la hauteur de lâcher
souhaitée. En posant l’hypothèse que le mouvement de la corde est quasi-statique
avant d’être lâchée, (9.1) devient :

F = KW. (9.2)
Afin de calculer F, nous devons renseigner le profil spatial de la corde avant son
lâcher. Puisque celle-ci est pincée ponctuellement en un nœud, son profil est trian-
gulaire avec pour sommet la valeur de hauteur de lâcher souhaitée au nœud pincé
(figure 14).
10-3
1

0.8
Hauteur (m)

0.6

0.4

0.2

0
10 20 30 40 50 60
Points discrets de la corde (noeuds)

Figure 14 – Profil de la corde avant d’être lâchée

36
Modèle non linéaire

Dans le cas non linéaire, la méthode des éléments finis conduit au système sui-
vant :
 
T
MẄ + CẆ + 1 + βW KW KW = F. (9.3)
Sous les mêmes hypothèses que dans le cas linéaire, le système devient :
 
F = 1 + βWT KW KW (9.4)
Le calcul de la force se réalise de la même manière que dans le cas linéaire, c’est à
dire en renseignant le profil spatial de la corde avant son lâcher dans le vecteur W.

9.2 Annexe 2 : Comparaison des différentes bibliothèques d’algèbre li-

néaire en C++
Afin de concilier performance et simplicité de transcription du code MATLAB en
C++, différentes bibliothèques dédiées à l’algèbre linéaire ont été testées. Les tests
ont eu lieu sur le modèle de corde linéaire discrétisée en 20 éléments. La charge du
CPU pour chacune des bibliothèques est donnée ci-dessous :
— It++ : 28%
— Eigen : 12%
— Modalys : 10%
— Coala : 50%
— BLAS/LAPACK : 14%
— module DSP de JUCE : 12%
— Armadillo : 13%
On constate alors que les libraires les plus performantes sont : Eigen, Modalys et le
module DSP de JUCE. Les deux derniers ne permettant pas une transcription aisée
du script MATLAB, Eigen est retenue dans ce projet.

37
Références
[1] Site du Master ATIAM. http://www.atiam.ircam.fr. Dernier accès :
03/07/2018.
[2] Site de l’IRCAM. https://www.ircam.fr. Dernier accès : 03/07/2018.
[3] Site de l’équipe S3AM. http://s3am.ircam.fr. Dernier accès : 03/07/2018.
[4] Site du logiciel Modalys. http://forumnet.ircam.fr/fr/produit/modalys/.
Dernier accès : 03/07/2018.
[5] J.L.R. d’Alembert. Recherches sur la courbe que forme une corde tendue mise
en vibration. Histoire de l’Académie royale des sciences et belles-lettres de Berlin
(HAB) pour l’année 1747. 1750.
[6] C. Cadoz. Synthèse sonore par simulation de mécanismes vibratoires. Applica-
tions aux sons musicaux. These, Institut national polytechnique de Grenoble,
Octobre 1979.
[7] C. Cadoz, A. Luciani, and J.-L. Florens. Responsive Input Devices and Sound
Synthesis by Stimulation of Instrumental Mechanisms : The Cordis System.
Computer Music Journal, 8(3) :60–73, Septembre 1984. Texte republié dans
l’ouvrage The music Machine (sous la direction de Curtis Roads), en 1989 :
pages 495-508.
[8] Site du logiciel CORDIS-ANIMA. http://www.acroe-ica.org/fr/
recherche/le-langage-de-modelisation-physique-cordis-anima. Der-
nier accès : 04/07/2018.
[9] J.-M. Adrien. The missing link : Modal synthesis. pages 269–298. MIT Press,
Cambridge, MA, USA, 1991.
[10] M. E. McIntyre, R. T. Schumacher, and J. Woodhouse. On the oscillations
of musical instruments. The Journal of the Acoustical Society of America,
74(5) :1325–1345, 1983.
[11] J. O. Smith III. Physical Audio Signal Processing : for Virtual Musical Instru-
ments and Digital Audio Effects. 2006.
[12] S. Bilbao. Numerical Sound Synthesis : Finite Difference Schemes and Simu-
lation in Musical Acoustics. Septembre 2009.
[13] K.J. Bathe. Finite Element Procedures. Prentice Hall, 1996.
[14] S. Bilbao and J. Smith. Energy-conserving finite difference schemes for nonlinear
strings. 91 :299–311, Mars 2005.
[15] J. Chabassier and P. Joly. Energy Preserving Schemes for Nonlinear Hamilto-
nian Systems of Wave Equations. Application to the Vibrating Piano String.
Research Report RR-7168, INRIA, 2010.
[16] B. Bank. Energy-based synthesis of tension modulation in strings. Octobre
2009.
[17] A. van der Schaft and D. Jeltsema. Port-Hamiltonian Systems Theory : An
Introductory Overview. Now Publishers, 2014.
[18] B.M. Maschke, A.J. Van Der Schaft, and P.C. Breedveld. An intrinsic hamil-
tonian formulation of network dynamics : non-standard poisson structures and
gyrators. Journal of the Franklin Institute, 329(5) :923 – 966, 1992.
[19] A. van der Schaft G. Golo, V. Talasila and B. Maschke. Hamiltonian discreti-
zation of boundary control systems. Automatica, 40(5) :757 – 771, 2004.

38
[20] R. Moulla, L. Lefèvre, and B. Maschke. Pseudo-spectral methods for the spatial
symplectic reduction of open systems of conservation laws. Journal of Compu-
tational Physics, 231(4) :1272 – 1292, 2012.
[21] A. Macchelli. Boundary energy shaping of linear distributed port-hamiltonian
systems. European Journal of Control, 19(6) :521 – 528, 2013.
[22] T. Hélie and D. Roze. Corde non linéaire amortie : formulation hamiltonienne
à ports, réduction d’ordre exacte et simulation à passivité garantie. In 13ème
Congrès Français d’Acoustique, Le Mans, France, Avril 2016.
[23] F. L. Cardoso-Ribeiro, D. Matignon, and L. Lefèvre. A structure-preserving par-
titioned finite element method for the 2d wave equation. IFAC-PapersOnLine,
51(3) :119 – 124, 2018. 6th IFAC Workshop on Lagrangian and Hamiltonian
Methods for Nonlinear Control LHMNC 2018.
[24] N. Lopes. Approche passive pour la modélisation, la simulation et l’étude d’un
banc de test robotisé pour les instruments de type cuivre. Theses, Université
Pierre et Marie Curie - Paris VI, 2016.
[25] A. Falaize. Modélisation, simulation, génération de code et correction de sys-
tèmes multi-physiques audios : approche par réseau de composants et formula-
tion Hamiltonienne à Ports. Theses, Université Pierre et Marie Curie - Paris
VI, July 2016.
[26] Manuel Goubert, Université Lille 1 - Onde progressive sinusoïdale dans
une corde ; relation de dispersion et vitesse de phase. http://res-nlp.
univ-lemans.fr/NLP_E_M02_G03_01/co/NLP_E_M02_G03_01.html. Dernier
accès : 12/07/2018.
[27] G. Carrier. On the nonlinear vibration problem of the elastic string. Quarterly
of Applied Mathematics, pages 157–165, 1945.
[28] S. S. Rao. The Finite Element Method in Engineering. Janvier 1982.
[29] M. Géradin and D. Rixen. Mechanical Vibrations : Theory and Application to
Structural Dynamics. Wiley, 1997.
[30] V. Manet. La Méthode des Éléments Finis : Vulgarisation des aspects mathé-
matiques, Illustration des capacités de la méthode. Lecture, August 2013.
[31] T. Hélie, A. Falaize, and N. Lopes. Systèmes Hamiltoniens à Ports avec approche
par composants pour la simulation à passivité garantie de problèmes conservatifs
et dissipatifs. In Colloque National en Calcul des Structures, volume 12, Giens,
France, May 2015.
[32] V. Duindam, A. Macchelli, S. Stramigioli, and H. Bruyninckx. Modeling and
Control of Complex Physical Systems : The Port-Hamiltonian Approach. Jan-
vier 2009.
[33] G. Golo. Interconnection structures in port-based modelling : Tools for analysis
and simulation. Janvier 2002.
[34] F. Cardoso Ribeiro. Port-hamiltonian modeling and control of a fluid-structure
system, 2016. Thèse de Doctorat, Université de Toulouse, Institut Supérieur de
l’Aéronautique et de l’Espace.
[35] T. Hélie and D. Roze. Sound synthesis of a nonlinear string using Volterra
series. Journal of Sound and Vibration, 314(0) :275–306, 2008.
[36] G. Kirchhoff. Vorlesungen über mathematische physik : Mechanik. 1877.

39
[37] Y. Wu and T. Shen. Finite element approximation of the nonlinear kirchhoff
string with boundary control. Juillet 2014.
[38] N. Lopes, T. Hélie, and A. Falaize. Explicit second-order accurate method
for the passive guaranteed simulation of port-hamiltonian systems. IFAC-
PapersOnLine, 48(13) :223 – 228, 2015. 5th IFAC Workshop on Lagrangian
and Hamiltonian Methods for Nonlinear Control LHMNC 2015.
[39] S. Aoues. Energy preserving discretization of port-Hamiltonian systems. Theses,
INSA de Lyon, Décembre 2014.
[40] Site de l’équipe AAS. https://www.applied-acoustics.com. Dernier accès :
23/07/2018.
[41] Site de l’équipe MODARTT. https://www.pianoteq.com. Dernier accès :
23/07/2018.
[42] Site de la société ARTURIA. https://arturia.com/. Dernier accès :
26/07/2018.
[43] Site du produit MATLAB. https://fr.mathworks.com/products/matlab.
html. Dernier accès : 24/07/2018.
[44] Site du framework JUCE. https://juce.com. Dernier accès : 26/07/2018.
[45] Site de la librairie Eigen. http://eigen.tuxfamily.org. Dernier accès :
24/07/2018.
[46] Transcriptions MATLAB/Eigen. https://eigen.tuxfamily.org/dox/
AsciiQuickReference.txt. Dernier accès : 24/07/2018.

40