Marc Judy Pierre

Alexis Colette

Keywords: systèmes Hamiltoniens à ports, acoustique, modèles à constantes localisées, conduit vocal. vi

Cette thèse traite de la modélisation et de la simulation de l'appareil vocal, et plus particulièrement du conduit vocal articulé. L'objectif principal est de décrire de manière concomittante la propagation des ondes acoustiques et l'écoulement de l'air en interaction avec les tissus et les muscles des voies aériennes. Cette thèse combine des hypothèses, des approches et des outils issus de la mécanique des fluides et de la thermodynamique pour produire un modèle non linéaire à constantes localisées. Ce modèle prend en compte la propagation acoustique, la présence de dérivations (pour la cavité nasale par exemple) et les mouvements dus aux phénomènes articulatoires. Ce travail est effectué dans le domaine temporel, dans le cadre du formalisme des systèmes Hamiltoniens à ports (sHp) qui garantit l'obtention de modèles passifs et énergétiquement bien posés. La construction du modèle à constantes localisées se déroule suivant plusieurs étapes : (a) partition du conduit vocal en tronçon élémentaires; (b) pour chaque tronçon, décomposition des champs de vitesse et de densité de masse sur des fonctions de forme associées à l'écoulement axial, à l'écoulement transverse et à la compression du fluide ; (c) projection des équations de conservation de la masse et de la quantité de mouvement, puis formulation Hamiltonienne à ports (Hp) macroscopique ; (d) interconnexion de ces modèles sHp à l'aide de leurs ports (avec de possibles dérivations) afin de construire le modèle non linéaire complet du conduit vocal. Nous démontrons, à l'aide de circuits électriques équivalents, que ce modèle apporte de nouvelles nouvelles interprétations sur les interactions fluides-structures (IFS) tout en préservant la modularité, l'interprétabilité et l'extensibilité des modèles de l'état de l'art. Pour obtenir une formulation acoustique et améliorer le conditionnement numérique, nous proposons un sHp équivalent dit «relevé »dont les variables sont exprimées, à l'aide d'un changement de variable, comme des fluctuations autour d'un état au repos. Le modèle complet du conduit vocal résulte de l'interconnexion de plusieurs tronçons et de modèles mécaniques des tissus incluant des parois mobiles, pilotées pour modifier la géométrie. Ces interconnexions multiples mettent en lumière la présence de relations algébriques (contraintes) que nous résolvons à l'aide d'une méthode basée sur les sHp différentiel-algébriques (sHp-DA) et les graphs dirigés. Nous démontrons, à l'aide de simulations à passivité garantie pour des coarticulations simples, l'impact de notre modélisation du mouvement des articulateurs sur le signal de sortie audible, tout en reproduisant les comportements attendus (bilan de puissance, résonances acoustiques, transitoires, convection de masse). Nous utilisons la même méthodologie de modélisation pour proposer de nouveaux modèles de jonction (de complexité différentes) entre trois tronçons ou tubes. Ces modèles peuvent être appliqués à l'acoustique des instruments à vent. Enfin, nous proposons une formulation Hamiltonienne à ports de fonctions indicatrices (ou de couleur) et des méthodes dites level-set pour la modélisation de sHp dont le domaine spatial dépend du temps. Cette approche nous permet de réinterpréter l'interconnexion et le couplage de système dynamiques dans le cas des interactions fluides-structures. v MOTS-CL ÉS: Interactions fluide-structure,

A B S T R A C T

This work addresses the physical modelling and simulation of the vocal apparatus with a focus on the articulated vocal tract. The objective is to model both the propagation of acoustic waves and the fluid flow interacting with the tissues along the upper airways. Getting back to the basic principles of fluid mechanics and of thermodynamics, this thesis combines a set of hypotheses, approaches and tools to produce a nonlinear lumped-parameter model that accounts for acoustic propagation, but also for side branching (e.g., nasal coupling) and for articulation (with hard and soft moving tissues). This work is carried out in the time-domain using the framework of port-Hamiltonian systems (pHs) that ensures the energy consistency and passivity of the models. The lumped-parameter model is built in several steps: (a) partition of the vocal tract into elementary tracts; (b) for each tract, decomposition of velocity and density fields on basis functions associated with the axial flow, transverse expansion flow, and fluid compression; (c) formulation of the projected mass and momentum conservation equations as a macroscopic pHs; (d) interconnection of these pHs through ports (possibly with side-branching) to build the nonlinear pHs model of the full vocal tract. Using electrical equivalent circuits, we show that this model provides new interpretations on Fluid-Structure Interactions (FSI) while preserving the modularity, physical interpretability and scalability of classical lumped-parameter models of the state-of-the-art.

To obtain an acoustical formulation an improve the numerical conditioning, we propose an equivalent (shifted) pHs, which variables are expressed in terms of fluctuations around a state at rest, using a variable change. The full vocal tract is obtained by interconnecting several tracts (using ports) and mechanical models of tissues on mobile walls, the latter being used to drive the geometry. This interconnection step involves algebraic relations that we solve with an assembly method based on directed graphs and differential-algebraic pHs. Power-balanced numerical simulations for simple co-articulations demonstrate the significant impact of the energy-consistent modelling of time-varying tracts on audible signals, while reproducing expected properties (power balance, acoustical resonances, transient effects, mass convection). As a corollary work based on the same hypotheses and approach, we propose lumped passive models of junction (of various complexity) between three tracts or pipes, that could be also used in musical acoustics for wind instruments. At last, we revisit the theoretical framework of the so-called indicator (or color) functions and level set methods for the time-space modelling of pHS with time-varying domains. This provides new insights on the interconnection and the coupling of infinite-dimensional FSI systems.

KEYWORDS: fluid-structure interactions, port-Hamiltonian systems vocal tract, acoustics, lumped-parameter modelling. Front cut of the vocal apparatus. . . . . . . . . . . . . . . . . . . . . . . Figure 2 Front cut of the larynx . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 3 A vocal tract approximated by concatenated cylinders. . . . . . . . . . . Figure 4 A linear mass-spring-damper (MKA) system. . . . . . . . . . . . . . . . Figure 5 A two-mass-spring-damper system where m 1 and m 2 are rigidly connected. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 6 A one-dimensional acoustical pipe with rigid and motionless walls. . . . Figure 7 Numerical simulation results for the mass-spring-damper system. . . . . Figure 8 Numerical results for different version of the discrete gradient (DG) method. Figure 9 Example of the shift procedure for a one-dimensional Hamiltonian function. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 10 Simplified sagittal cut of the vocal tract . . . . . . . . . . . . . . . . . . Figure 11 The lumped-parameter approach. . . . . . . . . . . . . . . . . . . . . . . Figure 12 Considered junction between K tracts through a common interface . . . Figure 13 2D cut of a tract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 14 Velocity field decomposition. . . . . . . . . . . . . . . . . . . . . . . . . . Figure 15 Equivalent electrical circuit of tract model . . . . . . . . . . . . . . . . . Figure 16 Equivalent electrical circuit of the tract model of Maeda . . . . . . . . . Figure 17 A full vocal tract shown as a network of interconnected tracts . . . . . . Figure 18 A basic two-tract example. . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 19 Oriented graph of the network in Fig. 17 . . . . . . . . . . . . . . . . . . Figure 20 Example: directed graph of the two-tract case . . . . . . . . . . . . . . . Figure 21 Example of a vocal tract with N t " 5. . . . . . . . . . . . . . . . . . . . . Figure 22 Directed graph of example 2 . . . . . . . . . . . . . . . . . . . . . . . . . Figure 23 Spring-Damper system . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 24 Equivalent electrical circuit of the shifted tract model augmented with simple resistors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 25 Two-tract system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 26 Height and total mass versus time at the dynamical phase. . . . . . . . . Figure 27 Local error on the powerbalance for the two-tract model. . . . . . . . . . Figure 28 Frequency responses of the two-tract model for both static phases . . . . Figure 29 Frequency response of a uniform, closed-open, 17-centimeter-long tube discretised with N t " 1, 22 tracts. . . . . . . . . . . . . . . . . . . . . Figure 30 Simulated geometry of the vocal tract. . . . . . . . . . . . . . . . . . . . Figure 31 Frequency response of the vocal tract for phases (1) and (3). . . . . . . Figure 32 Power flows of each tract during step (1) and (2). . . . . . . . . . . . . . Figure 33 Flow rate at the lips and its time derivative for a simple coarticulation. . Figure 34 Numerical results for the coarticulation a-to-i . . . . . . . . . . . . . . . Figure 35 Input enthalpy and output mass rate at both extermities of the vocal tract. Figure 36 A generic three port junction between 3 systems. . . . . . . . . . . . . . Figure 37 The five different junction models . . . . . . . . . . . . . . . . . . . . . .

C O N T E N T S

Considered geometry of the junction. . . . . . . . . . . . . . . . . . . . . Figure 39 Equivalent electrical electrical circuit of model (M1b) . . . . . . . . . . . Figure 40 A simple fluid-structure interaction problem. . . . . . . . . . . . . . . . Figure 41 The control domain Ω paved with subdomains ω i , i P 1, N ω . . . . . Figure 42 The color function c i and its gradient along the associated subdomain ω i Figure 43 Numerical results for the mass-spring-damper . . . . . . . . . . . . . . . Figure 44 Numerical results for the adimensioned mass-spring-damper system. . . . Figure 45 Velocity input signal for the second tract . . . . . . . . . . . . . . . . . . Figure 46 Enthalpy impulse used in the simulations . . . . . . . . . . . . . . . . . . Figure 47 A simple coarticulation: velocity signals applied to the walls. . . . . . . . Figure 48 Vowel This thesis is part of a collaboration between the S3AM 1 team of the STMS 2 laboratory housed by the IRCAM 3 in Paris, and the Sounds 4 team of the LMA 5 laboratory in Marseilles. Most of the work was carried out at the STMS laboratory.

the s3am team -The physical modelling of audio and acoustic systems is an essential topic of research for the S3AM team. Some of their recent research subjects include: the loudspeaker [START_REF] Falaize | Passive Modelling of the Electrodynamic Loudspeaker: From the Thiele-small Model to Nonlinear Port-hamiltonian Systems[END_REF], the Fender Rhode [START_REF] Falaize | Guaranteed-passive Simulation of an Electromechanical Piano: A Port-hamiltonian Approach[END_REF], the Ondes Martenot [START_REF] Najnudel | Analog circuits and Port-Hamiltonian realizability issues: a resolution method for simulations via equivalent components[END_REF], brass instruments [START_REF] Lopes | Energy Balanced Model of a Jet Interacting With a Brass Player's Lip[END_REF], analog audio circuits [Naj+21; MH19; FH13], the nonlinear string [START_REF] Hélie | Corde non linéaire amortie : formulation hamiltonienne à ports, réduction d'ordre exacte et simulation à passivité garantie[END_REF], the tom-tom drum [START_REF] Wijnand | Active control of the axisymmetric vibration modes of a tomtom drum[END_REF], and the vocal apparatus [START_REF] Hélie | Self-oscillations of a vocal apparatus: A porthamiltonian formulation[END_REF]. The broad range of applications is unified by the increasing usage of the port-Hamiltonian formalism [START_REF] Maschke | Port-Controlled Hamiltonian Systems: Modelling Origins and System Theoretic Properties[END_REF][START_REF] Van Der Schaft | Hamiltonian Formulation of Distributedparameter Systems with Boundary Energy Flow[END_REF], which allows for a systematic representation and energy consistent approach of physical modelling. It has become a common language between team members. The growing interest in port-Hamiltonian systems led to the PhD thesis of Antoine Falaize [START_REF] Falaize | Modélisation, simulation, génération de code et correction de systèmes multi-physiques audios: Approche par réseau de composants et formulation Hamiltonienne à Ports[END_REF], who created a Python library called PyPHS [Fal], dedicated to the modelling and simulation of port-Hamiltonian systems.

The S3AM team also works on Volterra Series [START_REF] Bouvier | Identification de systèmes non linéaires représentés en séries de Volterra : applications aux systèmes sonores[END_REF], finite-time control [Wij+18; Wij+19; WdR21], and differential geometry [START_REF] Bensoam | Geometric Numerical Methods with Lie Groups[END_REF] applied to physical modelling sound synthesis. The S3AM team was recently involved in three collaborative research projects: the two ANR projects Hamecmopsys 6 and finite4SOS 7 , and the European project INFIDHEM 8 .

the sounds team -The topics of research of the Sounds team are related to the audible frequency range: environment and noise pollution, auditory perception, and musical acoustics. They collaborate with major companies of their respective domains: EDF, Buffet Crampon, and Yamaha.

covid-19 -The thesis was carried out from October 2018 to December 2021. From March 2019 to the handing in of the manuscript, the COVID outbreak impacted our working conditions on several levels. Firstly, over the lockdown periods, most of the research was carried outside of the laboratory. Secondly, in August 2020, I caught the COVID-19 virus and was ill for three months, followed by three months of not being physically able to work full time. I am grateful to Sorbonne-Université for extending my contract by three months, which enabled me to carry my PhD project to term.

research project

The vocal apparatus is a complex multiphysical system capable of producing a great number of different sounds. These sounds are the result of complex fluid-structure interactions and acoustical phenomena. The main objectives of the collaboration between the S3AM team and the Sounds team, which this thesis is a part of, are:

1) to provide realistic physical models of the multiphysical components that form the vocal apparatus;

2) to produce power-balanced simulations of the whole system;

3) to recreate and control the known behaviours of both healthy and pathological voices.

The outcomes of this research project are targeted at the medical and artistic domains.

To ensure that the proposed nonlinear physical models guarantee the conservation of the energy, this research is carried out in the framework of port-Hamiltonian systems [START_REF] Maschke | Port-Controlled Hamiltonian Systems: Modelling Origins and System Theoretic Properties[END_REF][START_REF] Van Der Schaft | Port-Hamiltonian Systems Theory: An Introductory Overview[END_REF].

In this thesis, we consider objectives 1 and 2, and focused mainly on the vocal tract. This PhD project is the continuation of the work of Nicolas Lopes [START_REF] Lopes | Energy Balanced Model of a Jet Interacting With a Brass Player's Lip[END_REF], who, amongst other research topics, worked on a fluid-structure interaction model of the lips of a brass player.

outline of the thesis

The manuscript is organised as follows.

In chapter 1, we describe the subject of study and the physics involved in speech production. Then, we give an overview of the literature, describe our approach based on fluid mechanics and thermodynamics, and list our hypotheses with the associated starting physical equations.

In chapter 2, we introduce the port-Hamiltonian framework for the finite-dimensional and infinite-dimensional cases. We adopt a component-driven approach where passive physical systems are described as the interconnection of energy-storing dynamical components, passive (conservative or dissipative) memoryless components, and external ports that enable the system to interact with its environment. In a final section, we remind the reader of an existing port-Hamiltonian-compatible numerical method, then improve it. To shift from a multiphysical description of the vocal tract to a (nonlinear) acoustical one, we introduce a method that formulates a port-Hamiltonian system as fluctuations around a state at rest.

In chapter 3, we propose a new passive model of the vocal tract. To account for the energetic contribution of the moving walls, driven by the articulators, we adopt a two-dimensional fluid mechanical setting. To produce a simple model, we consider the class of lumped-parameter models, a popular approach for the physical modelling of the vocal tract. In this approach, the vocal tract is modelled as an interconnected network of elementary tracts of simple geometry. In chapter 4, we perform numerical experiments to examine the capacities and limitations of the proposed model and numerical method.

In chapter 5, we use the same hypotheses and modelling methodology as for the vocal tract to propose a new model of a three-port junction that is compatible with fluid mechanics and acoustics. This model can be used to describe branchings in the vocal apparatus (e.g. the nasal branching), but also to describe branchings in wind musical instruments (e.g. branched resonator, lateral holes).

Finally, in chapter 6, we introduce a new tool for the modelling of infinite-dimensional fluid-structure interactions problems. Based on the work of Diagne and Maschke [START_REF] Diagne | Port Hamiltonian Formulation of a System of Two Conservation Laws with a Moving Interface[END_REF], we propose a distribution-based formulation of the indicator function (also called color function). We then include the color functions in an infinite-dimensional port-Hamiltonian system to account for time-varying domains.

Part I M O D E L L I N G A P P R O A C H A N D T O O L S

D E S C R I P T I O N A N D P H Y S I C A L M O D E L L I N G O F T H E V O C A L A P PA R AT U S

In this chapter, we describe the research subject and define the scope of this thesis. In section 1.1, we describe the vocal apparatus and some of the fluid mechanical phenomena involved in the production of speech. Then, in section 1.2, we propose an overview of the state of the art in the physical modelling of the vocal apparatus. Finally, in section 1.3, we determine the scope of this thesis and list our hypotheses and associated equations.

the vocal apparatus

The term vocal apparatus denotes the organs involved in the production of vocal sounds. It is usually separated into three groups, indicated with solid colors on Fig. 1. The respiratory group is mainly comprised of the lungs, the diaphragm, and the trachea. It is the main pressure source of the vocal apparatus. During the inhaling phase, the diaphragm contracts, expanding the lungs and letting air in. During the exhaling phase, the diaphragm relaxes, thus compressing the lungs and ejecting air through the trachea towards the larynx. The phonatory group is mainly composed of the larynx (see Fig. 2), an organ made of cartilages held together and controlled by muscles. The vocal ligament and the thyroarytenoid muscle, covered by the lamina propria and the epithelium, form the vocal folds, the tension and adduction of which are controlled by a set of muscles (not detailed here). The constriction formed by the two vocal folds, opposing themselves to the airflow coming from the lungs, is called the glottis. The term phonation refers to the emergence of periodic oscillations of the vocal folds interacting with the glottal flow. This phenomena produces sound waves that propagates along the vocal tract.

Glottis

The articulatory group is comprised of the upper airways (e.g. the pharynx, the oral and nasal cavities), also called the vocal tract, and the articulators (e.g. the tong, the mandible, the lips, the velum). It is delimited upstream by the larynx, and downstream by the lips and nose. The geometry of the air cavities vary with time under the action of the articulators. The interaction between the air cavities and the airflow, and the walls of the vocal tract can result in the creation of a sound source, a phenomena at the heart of the production of consonants.

the production of consonants -Consonants are produced by the narrowing (constriction) or the complete closure of the vocal tract. The location of the constriction is called the place of articulation. A consonant is said to be voiced if phonation occurs during its production. The opposite case is said to be voiceless.

Consonants are classified by voicing, place of articulation and type of constriction (completely or partially closed). For instance, the [s] consonant, an alveolar fricative, is produced by forcing air through a narrow opening between the upper and lower teeth, which creates a noise source. During this process, the vocal folds are held apart and do not vibrate, making the /s/ a voiceless consonant. [b] is a bilabial voiced stop consonant, where the lips close completely and the vocal folds vibrate.

Remark 1 (Pinktrombone: an online articulated vocal tract): Pinktrombone is an online interactive vocal synthesizer where the user can modify the shape of the vocal tract by moving the articulators themselves (e.g. tongue, velum, lips). It was created by Neil Thapen and is available at https: // dood. al/ pinktrombone/ .

the fluid mechanics behind speech production -The phenomena at the heart of speech production are well described within the scope of compressible (and incompressible) fluid mechanics and fluid-structure interactions. In the larynx, the sustained oscillations of the vocal folds create sound waves that propagate along the vocal tract. This auto-oscillating phenomena is the consequence of the coupling between the glottal jet (fluid) and the vocal folds (structure).

In the case of vowels, the vocal tract mostly acts as a resonator, where the acoustic waves, coming from the larynx, are modified as they propagate through the irregular geometry of the upper airways.

For partially closed consonants, the constrictions in the vocal tract may generate aeroacoustical sources. The airflow coming from the lungs separates from the walls of the vocal tract, thus creating a jet. Downstream of the constriction, the jet progressively becomes unstable until it becomes fully turbulent [START_REF] Kenneth | Airflow and Turbulence Noise for Fricative and Stop Consonants: Static Considerations[END_REF], dissipating most of its energy into heat (viscosity effects) and a small fraction of it into sound waves, that propagate along the vocal tract.

For stop consonants, the complete closure of the vocal tract (e.g. at the lips, the back of the tongue) obstructs the main airflow. Upstream of the constriction, the pressure builds up until the air pressure and the force applied by the articulators balance one another. At the sudden opening of the constriction, a resulting pressure wave propagates along the vocal tract.

During trills (like the [r] consonant), the tongue (structure) and the main airflow (fluid) interact to generate low frequency periodic oscillations. This brief overview of the fluid mechanical phenomena involved in the production of speech shows how fluid-structure interactions are a key element of speech production.

a short literature overview in the physical modelling of the vocal apparatus

Measuring the acoustical, fluid mechanical or biomechanical behaviour of vocal apparatus requires invasive captors and costly medical imaging equipement. Therefore, scientists have relied on physical models to better understand how the voice works.

The existing physical models are devoted to the study of each group of the vocal apparatus. Here, we focus on the phonatory and articulating groups. For both of them, we find low-dimensional models, using simplified geometries and physics to produce lightweight models, and high-dimensional models, in which the equations of physics are solved on a dense mesh to accurately reproduce the complex shape of the airways. The former models are mostly dedicated to the (real-time) synthesis of speech, as they are relatively easy to compute. The latter are favored for their realism, which comes at the cost of heavy computations.

Physical Modelling of the Phonatory Group

The phonatory group is classically modelled as an incompressible jet (the glottal jet) interacting with a deformable structure (the vocal folds). Finite-element-based models are mostly used to account for the complex geometry and biomechanical behaviour of the vocal folds. Ref. [START_REF] Alipour | Mathematical Models and Numerical Schemes for the Simulation of Human Phonation[END_REF] proposes an extensive review of this type of models. In Ref. [de +03], the use of a three-dimensional setting enables the authors to account for a complex geometry. Considering a two-dimensional setting makes it possible to reduce the amount of computations needed [START_REF] Balázsová | Comparison of numerical methods for the solution of viscous incompressible and low Mach number compressible flow[END_REF]. To further reduce the dimensionality, some simple two-mass models of the vocal folds can be used, which enables a focus on the computations of the dynamics of the fluid [START_REF] Sváček | FE Numerical Simulation of Incompressible Airflow in the Glottal Channel Periodically Closed by Self-sustained Vocal Folds Vibration[END_REF]. High-dimensional approaches allows for the study of energy transfers between the vocal folds and the flow, crucial to the understanding of the phonation process. For instance, Thomson, Mongeau, and Frankel [START_REF] Scott | Aerodynamic Transfer of Energy to the Vocal Folds[END_REF] showed that the alternating diverging-converging shape of the glottis is a key element of the self-oscillating process.

Low-dimensional models reduce to a minimum the number of unknown variables only to capture the most important features. For instance, Ishizaka and Flanagan [START_REF] Ishizaka | Synthesis of Voiced Sounds From a Two-Mass Model of the Vocal Cords[END_REF] synthesised voiced sounds with a model of the larynx containing only two mass-spring-damper systems, and Titze [START_REF] I R Titze | The Human Vocal Cords: A Mathematical Model[END_REF] built a model of the larynx coupling 16 masses. Whereas these models only account for transverse movements (perpendicular to the main flow), Adachi and Yu [AY05] developed a model accounting for an axial movement. These simple models are still being improved today, like the model of Elie and Laprie [START_REF] Elie | Extension of the single-matrix formulation of the vocal tract: Consideration of bilateral channels and connection of self-oscillating models of the vocal folds with a glottal chink[END_REF] which includes the presence of a glottal chink (asymetrical closure of the vocal folds).

Physical Modelling of the Vocal Tract

In the same fashion as for the phonatory group, two different trends can be identified. Seeking a high degree of realism, high-dimensional models rely on finite element methods (FEMs) to accurately reproduce the acoustical characteristics of the vocal tract, caused by its geometric intricacies. Applying finite-difference time domain to the vocal tract [START_REF] Takemoto | Acoustic analysis of the vocal tract during vowel production by finite-difference time-domain method[END_REF] is fairly uncommon. The geometrical data are acquired by the means of resonant magnetic imaging [TMK10; Bad+02; Han+07; Arn+19] that allows, using image processing, for the measurement of the complex three-dimensional shape of the vocal tract, including the side cavities (e.g. the piri-form fossae [START_REF] Takemoto | Acoustic Interaction between the Right and Left Piriform Fossae in Generating Spectral Dips[END_REF]) and nasal cavities. The precision of high-dimensional models allow for the extensive study of acoustical properties for the vocal tract [Han+07; VH Š08].

These models may become very costly to compute. To reduce the dimensionality of the models, Arnela and Guasch [START_REF] Arnela | Two-dimensional vocal tracts with three-dimensional behavior in the numerical generation of vowels[END_REF] assumed that the cross sections of the vocal tract are of elliptical shape, which enabled them to consider a two-dimensional setting. Neglecting the bend of the vocal tract has a minor impact on the acoustical resonances (about 5%), as shown by Sondhi [START_REF] Mohan | Resonances of a bent vocal tract[END_REF], which allows for further simplifications.

Low-dimensional models rely on stronger geometrical simplifications to produce simple models compatible with real-time synthesis. A classical approach [START_REF] Kröger | Articulatory synthesis of speech and singing: state of the art and suggestions for future research[END_REF] is to divide the vocal tract into smaller cylindrical cavities (in which plane waves are assumed to propagate), as depicted on Fig. 3, called tracts. The resulting onedimensional models can realistically reproduce the acoustical behaviour of the vocal tract up to 3 kHz [VH Š08]. Two main approaches are distinguishable.

T 1 T 2 T ... T n T n+1 T ... T N
Figure 3: A vocal tract approximated by concatenated cylinders. In transmission-reflection models, the pressure field in each tract is decomposed in forward and backward travelling plane waves. In lumped-parameter models, the physical variables in each tract are considered to be homogeneous (piecewise constant discretisation).

The transmission-reflection models [KL62; Lil85; Str00] decompose the acoustic field within a tract as a superposition of forward and backward traveling waves. At each tract-tract interface, an incident wave is scattered in a reflected and transmitted wave, the magnitude of which is determined by the physical parameters and geometry of the junction. These models have been extended to account for conical tracts by Välimäki and Karjalainen [START_REF] Välimäki | Improving the Kelly-lochbaum Vocal Tract Model Using Conical Tube Sections and Fractional Delay Filtering Techniques[END_REF] and flared tracts by Mignot, Hélie, and Matignon [START_REF] Mignot | From a model of lossy flared pipes to a general framework for simulation of waveguides[END_REF]. Mullen, Howard, and Murphy [START_REF] Mullen | Waveguide physical modeling of vocal tract acoustics: flexible formant bandwidth control from increased model dimensionality[END_REF] pushed this approach to its limit, by drastically increasing the number of tracts and considering the two-dimensional propagation of acoustic waves. These models are inexpensive to compute, which explaining their popularity in the musical acoustic domain. However, as pointed out by Kröger and Birkholz [START_REF] Kröger | Articulatory synthesis of speech and singing: state of the art and suggestions for future research[END_REF], these models cannot handle the lengthening of the tracts that occurs in simple coarticulations (e.g. [a]-to-[i], [a]-to-[u]), a major obstacle for their use for voice synthesis.

The lumped-parameter models [Tit73; FIS75; Mae82; ST95; Sto95] rely on the same geometric assumptions, but consider that the physical quantities in a tract are homogeneous. It is then possible to build a low-dimensional physical model, the dynamics of which is described by a set of ordinary differential equations (ODEs). To obtain numerical simulations, the equations have to be approximated by a timediscretising scheme. Maeda [START_REF] Maeda | A digital simulation method of the vocal-tract system[END_REF] proposed such a model paired with a real-time implementation. This model was continuously improved, like in 2008 by Mokhtari, Takemoto, and Kitamura [START_REF] Mokhtari | Single-matrix formulation of a time domain acoustic model of the vocal tract with side branches[END_REF], and in 2016 by Elie and Laprie [START_REF] Elie | Extension of the single-matrix formulation of the vocal tract: Consideration of bilateral channels and connection of self-oscillating models of the vocal folds with a glottal chink[END_REF]. The popularity of these models is due to their modularity and scalability which allows for the addition of: source-noise models [START_REF] Birkholz | Construction And Control Of A Three-Dimensional Vocal Tract Model[END_REF] based on signal processing research [START_REF] Narayanan | Noise source models for fricative consonants[END_REF]; turbulence loss models [START_REF] Birkholz | Simulation of losses due to turbulence in the time-varying vocal system[END_REF]; tongue models [START_REF] Elie | Simulating alveolar trills using a two-mass model of the tongue tip[END_REF]. Paired with an adequate articulatory model, they allow for the realistic synthesis of humain voice.

All the aformentioned models of the vocal tract rely on acoustical descriptions, therefore neglecting the energy contribution of the moving wall. Therefore, the power balance of these models are only satisfied for statical configurations, when the vocal tract is motionless (e.g. sustained vowels). This inspired some research around the power-balanced modelling of the larynx [HS17; Mor+18], or the use of a fluid-structure interaction (FSI) framework to account for the movement of the walls [START_REF] Guasch | A stabilized finite element method for the mixed wave equation in an ALE framework with application to diphthong production[END_REF]. However, theses approaches remain quite rare.

a fluid mechanical, power-balanced and low-dimensional physical modelling approach

In this thesis, we are interested in establishing new FSI lumped-parameter passive models of the vocal apparatus. We rely on a compressible fluid mechanical description so as to account for acoustical and fluid mechanical phenomena.

Scope of This Thesis

We are interested in producing new physical models that:

1. are simple and low-dimensional;

2. account for the energetical contribution of the articulation and side branching;

3. describe the nonlinear fluid mechanical and acoustical phenomena;

4. are passive in the sense that they verify a physical power balance.

To meet these requirements, we make the following choices. We consider the class of lumped-parameter models. We focus on the study of the simple fluid mechanical effects, discarding the phenomena associated with the vorticity of the fluid (turbulences, noise sources). Indeed, vorticity and aeroacoustic sources are very complex topics that are beyond the scope of this thesis.

Even though we consider a fluid mechanical setting, we will be using descriptions that embed acoustical phenomena. We consider a two-dimensional spatial domain to account for the transverse motion of the fluid induced by the side branches and the movement of the articulators.

Finally, to guarantee the passivity of the models, we will be using the port-Hamiltonian formalism, a modular and multiphysical modelling framework, that we introduce in chapter 2.

Considered Hypotheses and Equations

We now state the hypotheses of compressible fluid mechanics (H1-2) and thermodynamics (H3-5) used in this thesis.

Regarding the flows and fluids, we assume that:

(H1) the fluid is inviscid (no viscosity);

(H2) the flow is irrotational.

These hypotheses enable us to describe the acoustic and mass convection phenomena, and to neglect the presence of vorticity and turbulences.

Regarding the thermodynamical behaviour, we assume that:

(H3) the fluid (air) fullfils the perfect gas assumption;

(H4) the fluid undergoes an adiabatic and isentropic transformation (reversible process);

(H5) the heat capacity at constant volume C V is assumed to be constant w.r.t. the temperature within the considered range of fluctuations.

Under (H4), the heat sources are neglected. Note that this set of hypotheses does not assume any linearising assumptions so that it deals with nonlinear physics.

These hypotheses are linked to a set of conservation equations and thermodynamical laws. For a concise review on the mechanics of continua, we refer to the first chapter of Ref. [START_REF] Filippi | Acoustics : basic physics, theory, and methods[END_REF].

Consider a fluid contained in a virtual control volume Ω, equipped with a position vector r " px, y, zq P Ω. Given our hypotheses and assuming a Eulerian representation, a fluid is described by: its velocity field vpr, tq; its volumetric mass density ρpr, tq, also called mass density; its static pressure ppr, tq; its specific internal energy , expressed as per unit of mass; its temperature T ; its initial state characterised by the set of constant physical quantities tρ 0 , T 0 , P 0 u, where ρ 0 and P 0 are the rest value of the mass density and pressure of the atmosphere.

For concision, the space-time dependence of these physical variables is omitted in the rest of the thesis, except when required to avoid ambiguity.

The dynamics of the mass density (without internal sources) is given by the conservation of mass equation [START_REF] Leal | Advanced Transport Phenomena: Fluid Mechanics and Convective Transport Processes[END_REF] Bρ Bt " ´divpρvq.

Under (H1-3), the dynamics of the velocity field is given by the conservation of momentum equation

Bρv

Bt

`divpρv b vq `gradppq " 0.

The time variation of the energy of the system is given by the conservation of energy equation

ρ B Bt ˆ `1 2 v ¨v˙" ´divppvq,
where is described by the following laws and relations:

(Joule's law) d " C V dT , ( H3,5)
(perfect gas law, specific form) T " p {ρr 0 , (H3)

(adiabatic law) p " P 0 p ρ {ρ 0 q γ , ( H4
)
where C V is the heat capacity at constant volume, r 0 is the specific universal gas constant, γ " Cp {C V is the heat capacity ratio (" 1.4 for an ideal gas) with C p the heat capacity at constant pressure, and ρ 0 is the volumetric density at rest.

Finally, under (H3-5), the pressure is linked to the specific internal energy by the following relation [vM02] p " ρ 2 B pρq Bρ .

summary and perspectives

In this chapter, we briefly described the research subject and the physical phenomena involved in the production of speech. After proposing an overview of the state of the art, we stated our approach and the fluid mechanical phenomena we are tackling in this thesis. We ended the chapter by giving a compatible set of hypotheses and their associated equations and laws.

P O RT -H A M I LT O N I A N S Y S T E M S

The port-Hamiltonian (pH) framework extends the classical Hamiltonian mechanics to include dissipation phenomena and external interactions. Introduced by Maschke and van der Schaft [START_REF] Maschke | Port-Controlled Hamiltonian Systems: Modelling Origins and System Theoretic Properties[END_REF], port-Hamiltonian sytems (pHs) model physical systems using energy as a common denominator, which makes them particularily suited for multi-physical problems.

The pH formalism relies on a systematic geometric structure that guarantees the passivity of the model. This way of writing dynamical systems also brings physical insights, allows for comparisons between models, and physical analogies. The growing interest towards the pH formalism is due to its energetic guarantees, modularity and relevance in physical modelling and control problems. This formalism is used in electrical engineering [MSM04; MH19; Fal16; Naj+18], fluid mechanics [Car16; Lop16], solid mechanics [START_REF] Brugnoli | Port-Hamiltonian Formulation and Symplectic Discretization of Plate Models Part I: Mindlin Model for Thick Plates[END_REF], thermodynamics [START_REF] Tobias | An Object-Oriented Library for Heat Transfer Modelling and Simulation in Open Cell Foams[END_REF] acoustics [START_REF] Hélie | Corde non linéaire amortie : formulation hamiltonienne à ports, réduction d'ordre exacte et simulation à passivité garantie[END_REF][START_REF] Trenchant | Discrétisation et commande de systèmes vibro-acoustiques, une approche hamiltonienne à ports[END_REF] and chemistry [START_REF] Ramirez | Irreversible Port-hamiltonian Systems: A General Formulation of Irreversible Processes with Application to the CSTR[END_REF]. Applications range from the modelling of old and decaying electroacoustic instruments [START_REF] Najnudel | Analog circuits and Port-Hamiltonian realizability issues: a resolution method for simulations via equivalent components[END_REF] and district heating networks [START_REF] Hauschild | Port-Hamiltonian modeling of district heating networks[END_REF] to the modeling and control of plasma fluxes in fusion-based nuclear reactors [START_REF] Minh | Port-hamiltonian Formulation for Systems of Conservation Laws: Application to Plasma Dynamics in Tokamak Reactors[END_REF]. For an extensive overview on pH-based modelling, we refer to Ref. [START_REF] Rashad | Twenty years of distributed port-Hamiltonian systems: a literature review[END_REF].

In section 2.1, we remind the reader of the pH formalism in the finite dimensional setting. In section 2.2, we introduce the infinite-dimensional pH formalism for distributed parameter systems. Finally in section 2.3, we consider the time discretisation aspects. We remind the reader of a pH compatible numerical scheme, improve it, and then introduce a new method to formulate shifted pHs.

finite-dimensional port-hamiltonian systems

For the need of this thesis, we only introduce pHs as a modelling tool. For an extensive take on the subject, we refer to Refs. [START_REF]Port-Hamiltonian systems : an introductory survey[END_REF][START_REF] Van Der Schaft | Port-Hamiltonian Systems Theory: An Introductory Overview[END_REF].

We consider a component-based approach where a physical system is seen as a collection of basic components interconnected by a conservative interconnection, analogous to an electrical circuit made of discrete components (resistors, capacitors, inductors) connected by a printed circuit board (PCB).

Component-Based Modelling

The energy variation of an open physical system satisfies the power balance

dE dt `Pdiss `Pext " 0 (2.1)
where

• E is the energy stored by the system;

• P diss ě 0 is the dissipated power;

• P ext is the power given by the system to the exterior (receiver convention).

Physical systems verifying this equality are called passive [START_REF] Van Der Schaft | L2-gain and Passivity Techniques in Nonlinear Control[END_REF]Ch. 4]. Each phenomenon (energy storing, dissipation and sources) is described by its own class of component. energy storing dynamical components -They are described by n x energy variables x n gathered in the energy variable vector x " rx 1 , x 2 , . . . , x nx s (also called the state vector ). The energy stored by the system is given by the Hamiltonian Hpxq , a storage function bounded from below [vJ14]. Under these definitions and the chain rule, the first term of Eq. (2.1) now reads

dE dt " dHpxq dt " ∇ x Hpxq dx dt " e f
where we introduce the efforts e " ∇ x Hpxq and flows f " dx dt of the pHs.

memoryless algebraic components -They describe instantaneous transfers of power. The vector of n w memoryless variables w " rw 1 , . . . , w nw s is paired with a vector of algebraic laws zpwq " rz 1 pwq, . . . , z nw pwqs so the dissipated power P diss reads

P diss " zpwq w $ &
% " 0 for conservative memoryless components, ě 0 for dissipative memoryless components.

Note that this structure allows for the modelling of irreversible processes resulting in heat transfers [START_REF] Najnudel | Power-balanced Modeling of Nonlinear Coils and Transformers for Audio Circuits[END_REF]. Some algebraic laws of variables w (extracted from w) and laws zpwq (extracted from zpwq) can be written as a matrix times vector operation such as zpwq " Wpxqw where: Wpxq " ´Wpxq in the case of conservative algebraic components (transformers, gyrators), or where Wpxq is non-negative for linear dissipations (linear resistor, dashpot). external ports -At last, the system may interact with its environnement through n p external ports (that can be control ports). The power P ext given by the system and received by the exterior to the system (receiver convention) is the product of the input vector u " ru 1 , . . . , u np s and output vector y " ry 1 , . . . , y np s , that is

P ext " y u.
(2.2)

Under the above definitions, the full power balance (2.1) reads dE dt `Pdiss `Pext " 0 ùñ ∇ x Hpxq dx dt `zpwq w `y u " 0 .

interconnecting components -Components are connected through a conservative interconnection, called a Dirac structure1 , that guarantees the passivity of the system. In this thesis, we consider the class of pHs that can be written in the following differential-algebraic representation where a skew-symmetric matrix Spxq " ´Spxq links the flow vector F " rf , w , y s and effort vector E " re , zpwq , u s such as

¨f w y ‹ ‹ ' F " Spxq ¨e zpwq u ‹ ‹ ' E (2.3)
where S is a square matrix of size dim E " dim F. Each non-null (pair of) entry of S defines a routing of power between two components.

Here, the skew-symmetric property of S ensures the power-balance as 0 " E SpxqE (2.3) " E F. The passivity is guaranteed by the power-balance, by the lower bound on Hpxq, and by P diss ě 0.

The interconnection matrix is divided according to the dimensions of the storage, algebraic and external port components vectors such as

S " » - - - J xx G w G p ´G w J ww G wp ´G p ´G wp J pp fi ffi ffi fl (2.4)
The interconnection matrix is said to be canonical if its entries are independent of the energy variables x and the physical parameters of the pHs. In this thesis, the power exchanges between the components are then encoded by ˘1 (elementary connection). For the general case, please refer to Ref. [START_REF] Van Der Schaft | Port-Hamiltonian Systems on Graphs[END_REF].

Remark 2 (Classical Differential Form of PHS and Integrability Condition): When the algebraic part of Eq. (2.3) can be solved with respect to w, the system a can be written as an input-state-output pHs [START_REF] Maschke | Port-Controlled Hamiltonian Systems: Modelling Origins and System Theoretic Properties[END_REF] $ & % dx dt " pJpxq ´Rpxqq ∇ x Hpxq `Gp u ´y " G p pxq e `Jpp pxqu (2.5) where matrices J " ´JT and R " R ě 0 are computed from matrix S and laws zpwq. This differential form is associated to a geometric Dirac structure if J (with its depedency w.r.t. x) satisfies an integrability condition (not detailed here, see Ref. [vJ14, p. 48]). Otherwise, the system is called a pseudo-Hamiltonian system [START_REF] Alejandro | Port-Hamiltonian Modelling of Fluid-structure Interactions in a Longitudinal Domain[END_REF][START_REF] Van Der Schaft | Port-Hamiltonian Systems Theory: An Introductory Overview[END_REF]. The integrability condition on J involves the matrices J xx , J ww , and G wp (see Eq. (2.4)) and the conservative part of the memoryless laws zpwq in Eq. (2.3).

a The output ŷ " ´y is used here to obtain a generator convention: ŷ u is the power supplied to the system. example: the mass spring dashpot system -Consider the linear mass-springdashpot system depicted on Fig. 4, where a mass, a spring and a damper are rigidly connected. The linear spring, of stiffness k, is described by its elongation x k and its stored energy 1 2 kx 2 k . The mass m 1 is described by its velocity v 1 and kinetic energy 1 2 mv 2 1 . The dashpot, of elongation rate v d (velocity), is connected to the mass and produces a reaction force equal to av. At the action point (denoted in red on the figure), the mass is submitted to the external force F ext , the power-conjugated output of which is the velocity v ext of the action point. The rigid connection is expressed as 9

x k " v d " v ext " v 1 where 9

x k is the elongation rate of the spring, defined as the difference of velocity between both its extremities. As the point of action is attached to the mass, the external velocity and the one of the mass are equal. As the damper and spring are attached to the mass on one side, and to a rigid frame on the other, their elongation rate are equal to the velocity of the mass.

We now give the port-Hamiltonian (pH) formulation of this simple dynamical system. First, we model the spring, mass, damper and action points as individual components. Then, we connect them together. We define the state vector x " rπ 1 , x k s where π " m 1 v 1 is the momentum of the mass. The Hamiltonian is the sum of both the kinetic and the potential energies expressed in terms of the energy variables x, such as

Hpx " rπ 1 , x k s q " 1 2m 1 π 2 1 `1 2 kx 2 k leading to the efforts e " ∇ x Hpxq " ˜π1{m 1 " v 1 kx k " F k ¸,
where F k is the reaction force of the spring.

The damper is described by the memoryless component of variable w " v d and linear algebraic law zpwq " aw regenerating the reaction force of the damper F d . The interaction at the action point is described as an external port of input u " ´Fext (reaction force received by the exterior) and output v ext .

Finally, the algebraic-differential formulation of the system reads acceleration of the mass elongation-rate of the spring elongation-rate of the dashpot velocity of the action point where the first row is the sum of forces acting on the mass, and the first column encodes the rigid connection between the mass, the damper and the spring. The matrix is divided according to Eq. (2.4). Note that with this choice of energy variables, the interconnection matrix is canonical.

¨9 π 1 9 x k w v ext ‹ ‹ ‹ ‹ ' " » - - - - - 0 ´1 ´1 ´1 1 0 0 0 1 0 0 0 1 0 0 0 fi ffi ffi ffi ffi fl ¨π1 {m 1 x k k aw ´Fext ‹ ‹ ‹ ‹ '

Interconnecting pHs

The class of pH systems are stable by interconnection, meaning that the conservative interconnection of two pHs remains a pHs (see Ref. [vJ14, Sec. 6.2, p. 72]). Consider two pHs of energy variables x 1 and x 2 with inputs u 1 and u 2 , outputs y 1 and y 2 and the subsets ūi"t1,2u Ď u i"t1,2u and ȳi"t1,2u Ď y i"t1,2u that will be connected together. We set N " dim u 1 " dim u 2 (" dim y 2 " dim y 2 ).

There are two ways of interconnecting pHs.

In the port-based connection (compatible connection), the input/ouput relations can be written as a conservative interconnection

« ū1 ū2 ff " αpx 1 , x 2 q « 0 A ´A 0 ff « ȳ1 ȳ2 ff ,
where α is a coupling factor, that may be function of the states and physical parameters.

As the interconnection matrix is skew-symmetric, the interconnection is conservative. The port-based connection can be used only if the external port configuration of both pHs is of compatible causality, meaning that the output of system 1 can be directly connected to the input of system 2, and vice versa. For instance, a mass and a spring are compatible as their causality are complementary (see the above linear mass-spring-damper example). Using a compatible connection restores a pHs, the formulation of which is given in Eq. (2.3).

In the interconnection by energy [vJ14] (incompatible connection), the causality of two pHs is identical leading to the following formulation

0 " A « ȳ1 ȳ2 ff
where A is of size N ˆp2N q, (2.7) defining a set of N constraints, as the dynamics of both systems needs to be augmented by a set of Lagrange multipliers λ in order for the constraints to be continuously satisfied.

DA-pHs

A differential-algebraic port-Hamiltonian system (DA-pHs) [van13] is a pHs of form (2.3) augmented by N λ Lagrange multipliers, gathered in the vector λ, and N λ algebraic constraint equations expressed as a function of the efforts e, such as

¨f w 0 y ‹ ‹ ‹ ‹ ' " » - - - - - - - J G w C λ G p ´G w J ww 0 G wp ´C λ 0 0 0 ´G p ´G wp 0 J pp fi ffi ffi ffi ffi ffi ffi fl ¨e zpwq λ u ‹ ‹ ‹ ‹ ' (2.8)
where C λ is the constraint matrix (of size n x ˆNλ ) encoding Eq. (2.7) and 0 are matrices, of appropriate size, filled with zero entries.

Note that the descriptor-pHs [START_REF] Beattie | Port-Hamiltonian descriptor systems[END_REF] framework also accounts for constraints and use of standard DAE theory [KM07; KM06] and is close to the Differentialalgebraic port-Hamiltonian (DA-pH) form used in this thesis.

We now use the mass-spring-damper as an example, connecting a second mass rigidly. example: the two-mass spring damper system -Consider the pH formulation comprised of the mass-spring-damper, shown on Eq. (2.6), and a second mass m 2 , of velocity v 2 and kinetic energy 1 2 m 2 v 2 2 (see Fig. 5). Both masses are rigidly connected, so that v 1 " v 2 . The state vector and Hamiltonian read Hpx " rπ 1 , π 2 , x k s q " π 2 1{2m1 `π2 2{2m2 `1{2kx 2 k , leading to the the effort vector e " rv 1 , v 2 , F k s .

The DA-pHs formulation reads

rigid connection ¨9 π 1 9 π 2 9 x k w 0 v ext ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' " » - - - - - - - - - - - ´1 ´1 ´1 1 ´1 1 1 1 ´1 1 fi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi fl ¨π1 {m 1 " v 1 π 2 {m 2 " v 2 x k k " F k aw " F d λ F ext ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' constraint force
where ´C λ " r1, ´1, 0s encodes the constraints e π 1 ´eπ 2 " v 1 ´v2 " 0 and where the dynamics of both m 1 and m 2 are altered by the Lagrange multipliers. When the system is linear, it is usual to compute an equivalent mass in order to eliminate the constraint and get a classical unconstrained pHs.

To do so, one can use the equivalent component method [START_REF] Najnudel | Analog circuits and Port-Hamiltonian realizability issues: a resolution method for simulations via equivalent components[END_REF] or, in case of multiple constraints, its generalisation that projects the constrained pHs in a constrained subspace [START_REF] Maschke | Port-Controlled Hamiltonian Systems: Modelling Origins and System Theoretic Properties[END_REF][START_REF] Van Der Schaft | Port-Hamiltonian Differential-Algebraic Systems[END_REF][START_REF] Wu | Port Hamiltonian System in Descriptor Form for Balanced Reduction: Application to a Nanotweezer[END_REF][START_REF] Cardoso | Port-hamiltonian Modeling and Control of a Fluid-structure System : Application to Sloshing Phenomena in a Moving Container Coupled to a Flexible Structure[END_REF]. In the nonlinear case, it casts a fully implicit constrained pHs into a semi-explicit pHs. The latter method is used in chapter 3.

The pH formalism is a powerful modelling tool that systematically ensures the passivity of our models. We now introduce the infinite dimensional pHs, an extension of the pH formalism for distributed parameter system.

infinite-dimensional port-hamiltonian systems

The finite-dimensional pH framework has been extended to the infinite-dimensional case [vM02; Vil07; JZ12] to account for physical quantities that depend on space and time, the dynamics of which is governed by partial differential equation (PDE). Again, refer to Ref. [START_REF] Van Der Schaft | Port-Hamiltonian Systems Theory: An Introductory Overview[END_REF]Ch. 14] for an extensive take on the subject.

In this work, we only consider conservative equations. For the modelling of dissipative phenomena, we refer to [START_REF] Matignon | A Class of Damping Models Preserving Eigenspaces for Linear Conservative Port-hamiltonian Systems[END_REF]. For irreversible processes, see Irreversible port-Hamiltonian systems in Refs. [START_REF] Ramirez | Irreversible Port-hamiltonian Systems: A General Formulation of Irreversible Processes with Application to the CSTR[END_REF][START_REF] Alejandro | Port-Hamiltonian Modelling of Fluid-structure Interactions in a Longitudinal Domain[END_REF]. Consider a domain Ω with boundary BΩ on which we define a state vector X P L 2 pΩ, R n q, which depends on space through the position vector r P Ω. The Hamiltonian function is replaced by a smooth functional on Ω HpX q " ż Ω HpX q dV, where HpX q is an energy density function.

The Hamiltonian gradient, producing the efforts e, is replaced by the variational derivative δHpX q {δX of HpX q defined uniquely2 by HpX ` δX q " HpX q ` ż Ω δHpX q δX δX dΩ `Op 2 q , @ P R, with δX being a smooth variations of X .

A state X pr, tq, with Hamiltonian HpX q, has a conservative dynamics governed by BX Bt " J δHpX q δX , which also written as f " J e,

where J denotes a skew-adjoint differential operator (called linear Hamiltonian operator ) that satisfies the Jacobi identity and skew-symmetry property (see definition 7.1 in [START_REF] Olver | Applications of Lie Groups to Differential Equations[END_REF]p. 435]) and where f " BX Bt and e " δHpX q δX . At the boundary, the system interacts with the exterior through its boundary variables f B and e B , respectively defined as the trace of B t X and δ X H on BΩ. For a well-posedness definition regarding the input/output configuration, we refer to Refs. [START_REF] Le Gorrec | Dirac structures and Boundary Control Systems associated with Skew-Symmetric Differential Operators[END_REF][START_REF] Jacob | Linear Port-Hamiltonian Systems on Infinitedimensional Spaces[END_REF]. The power balance states that the variation of stored energy equals the power exchanged at the boundary

dHpX q dt " f B e B .
example: the webster horn equation -As an example, we take the case of the Webster horn equation, adapted from Ref. [START_REF] Matignon | A Class of Damping Models Preserving Eigenspaces for Linear Conservative Port-hamiltonian Systems[END_REF] to obtain a canonical formulation3 . e µ (0, t)• e v (0, t)

e µ (L 0 , t) • e v (L 0 , t) Ω L0 x = 0 0 S(x, t)
Figure 6: A one-dimensional acoustical pipe with rigid and motionless walls.

We consider the case of an axisymmetric acoustic pipe of length L, the cross section Spxq of which depends on the spatial coordinate x, as depicted on Fig. 6. The spatial domain is Ω " r0, Ls. We assume that the fluid obeys the hypothesis of linear acoustics, the dynamics of which is governed by the linearised conservation equations of mass and momentum

Bv Bt " ´grad ˆp ρ 0 ˙and Bρ Bt " ´1 Spxq divpSpxqρ 0 vq,
where v " vpx, tq is the acoustic velocity, ρ " ρpx, tq and p " ppx, tq are the acoustic mass density, ρ 0 and P 0 are the rest state mass density and pressure, and where divp ¨q and gradp ¨q simplify to B Bx in this one-dimensional case.

The energy E stored by the fluid is

E " ż Ω ˆ1 2 ρ 0 v ¨v `ρ pρq ˙Spxq dx
where pρq is the specific internal energy of the fluid. The pressure p is defined by p " ´ρ2 B pρq Bρ « c 2 0 ρ with c 0 " a γP 0{ρ 0 being the speed of sound for a perfect fluid for an isentropic process.

We now give the pH formulation. Choosing the linear mass density µ " µpx, tq " ρSpxq and the mean axial velocity of the fluid v as energy variables yields the following Hamiltonian HpX " rv, µs q " where we choose f B " r´e v p0, tq, e v pL, tqs and e B " re µ p0, tq, e µ pL, tqs as boundary variables.

ż Ω ˆ1 2 ρ 0 v ¨vSpxq `1 2 γP 0 Spxqρ 2 0 µ 2 ˙dx
Such a system can be spatially discretised using a pH compliant discretisation scheme, like the recent partitioned finite element method (PFEM) introduced by Cardoso Ribeiro, Matignon, and Lefèvre [START_REF] Luiz | A Partitioned Finite Element Method for power-preserving discretization of open systems of conservation laws[END_REF], in order to obtain a finite dimensional pHs. To perform numerical simulations of a finite dimensional pHs, the set of ODE then needs to be approximated in the discrete-time domain by numerical scheme.

numerical method and conditioning

The Discrete Gradient: A Passive Numerical Method for pHs

To preserve the passivity, ensured by the pH framework, in the discrete-time domain, we use a tailored numerical method that preserves the passivity of the simulations. While the main concern of classical time integrators is the time derivative operator, the discrete gradient (DG) method [START_REF] Itoh | Hamiltonian-conserving Discrete Canonical Equations Based on Variational Difference Quotients[END_REF][START_REF] Lopes | Approche passive pour la modélisation, la simulation et l'étude d'un banc de test robotisé pour les instruments de type cuivre[END_REF] where F s " δt ´1 is the sampling frequency.

For linear pHs, equipped with a quadratic Hamiltonian, the system is explicit. The Hamiltonian is then of form Hpxq " p 1 {2qx Qx and the efforts are e " ∇ x Hpxq " Qx. The corresponding discrete gradient simply reads ∇ d

x Hpx k , δx k q " Q px k `δx k{2q. In the nonlinear case, the value of δx k is found by solving the implicit problem with an iterative solver, such as the Newton-Raphson method [BV04, Sec. 9.5]. However, the class of quadratic pHs can be manipulated, using a quadratisation method [Lop16, Sec. 3.3.1, p. 107], to render the problem explicit.

Forma [START_REF] Forma | Étude d 'un système non linéaire : la voix[END_REF] showed, through a series of numerical experiments on a simple model of the larynx [HS17; HSW19], the advantages of a passive numerical method over classical Runge-Kutta methods (e.g. ode23s, ode15s and ode23tb implemented in Matlab).

The DG method is implemented in PyPHS, an open-source Python library that we introduce hereafter.

PyPHS: An Open Source Library for Assembling and Simulating pHs

PyPHS [Fal] is an open-source Python library dedicated to the manipulation and simulation of pHs. It is build and maintained by Antoine Falaize [START_REF] Falaize | Modélisation, simulation, génération de code et correction de systèmes multi-physiques audios: Approche par réseau de composants et formulation Hamiltonienne à Ports[END_REF] 5 .

A pHs is implemented by defining the energy variables and the Hamiltonian as symbolic variables. Then, the user specifies the interconnection matrices, the physical and numerical parameters from which PyPHS generates a C++ or Python implementation of the simulation code.

On Fig. 7, we present some numerical results for simple experiment with the MKA example. The full code is listed in App. A.1. Here: m " 0.5 kg, k " 4441 N ¨m´1 so that f " 15 Hz; a " 5 kg ¨s´1 . On the top, the time evolution of the energy variables. In the middle the time evolution of the input F ext . On the bottom left, the phase space, and on the bottom right the power balance.

The power balance is verified as the power supplied to the system equals the sum of the time variation of stored energy and dissipated power (see the bottom right plot on Fig. 7).

However, the discrete gradient may be prone to cancelation errors. We demonstrate this phenomenon with a simple nonlinear Hamiltonian and show how to improve the numerical conditioning of the method.

Improved Discrete Gradient

Cancelation errors may arise at a given current state x ‹ k around which the increments δx k are small. In this situation, both the substraction and division operations in the discrete gradient are very sensible to cancelation errors. We remind the reader of this discrete operator (for separable Hamiltonian) and highlight both problematic operations in the following equation

r∇ x Hpxqs n « " ∇ d x Hpx, δxq ı n " $ ' ' ' & ' ' ' % h n px n `δx n q ´hn px n q δx n if |δx n | ą DG , h 1 n px n q otherwise,
where we remind the reader that is the numerical precision of the chosen floating point representation.

On Fig. 8, we plot the absolute value of the difference between a given discrete gradient operator and the value of the continuous gradient at the state x ‹ k . When δx k ď DG , both the numerator and the denominator go to zero, leading to round off errors in the substracting and dividing operator resulting in an increasingly large error as δx k Ñ 0 (see the blue curve for δx ă ? on Fig. 8a). The desired behaviour is for the discrete gradient to keep converging to zero as δx k Ñ 0.

To fix this numerical problem, we adapt DG , enabling the discrete gradient to switch from the finite difference approximation to the alternate expression before the cancellation errors start increasing. In the case of Fig. 8,setting 6 DG " ? enables the switch to the other approximation to happen sooner, leading to a monotonic variation of the cancellation errors (see the orange curve on Fig. 8a).

However, the discontinuity that was present for δx k " now is placed at δx k " ? , which can be a problem if a dynamical system exhibits small variations of δx k of the same magnitude as DG . This behaviour can be misleading for numerical iterative solvers, as this discontinuity may be often triggered.

To solve this undesired behaviour, we extend the Taylor developpement of the continuous approximation, allowing for a smoother transition and a better approximation of the gradient of the Hamiltonian (dotted green curve on Fig. 8b). The improved discrete gradient (IDG) then reads

r∇ x Hpxqs n « " ∇ d x Hpx, δxq ı n " $ & % hnpxn`δxnq´hnpxnq δxn if |δx n | ą DG , h 1 n px n q `h2 n px n qδx n otherwise.
Note that for non-separable Hamiltonians, a similar modification of the discrete gradient can be proposed in Ref. [START_REF] Falaize | Modélisation, simulation, génération de code et correction de systèmes multi-physiques audios: Approche par réseau de composants et formulation Hamiltonienne à Ports[END_REF][START_REF] Hairer | Geometric numerical integration: structure-preserving algorithms for ordinary differential equations[END_REF] by using the Hessian of H in place of h 2 n (diagonal case).

10 -13 10 -10 10 -7 10 -4 10 -1 δx 10 -16 10 -13 10 -10 10 -7

10 -4

10 -1 10 2 ||∇ d x H(x k , δx k ) -∇ x H(x k )|| √ ε √ ε ε ε Below machine's precision

Cancelation errors

Discrete gradient, Eq. (1.11)

DG = √ DG =
(a) Impact of the value of DG on the DG.

10 -13 10 -10 10 -7 10 -4 10 -1 δx 10 -16 10 -13 10 -10 10 -7

10 -4 Remark 3 (Adimensionned models): In the case of an adimensioned nonlinear pHs (see App. A.2), DG should be mutipled by the inverse of the scaling factor.

10 -1 10 2 ||∇ d x H(x k , δx k ) -∇ x H(x k )|| √ ε √ ε ε ε Below machine's precision

Shifting pHs Around a State at Rest

The state at rest x ‹ of a closed dynamical system dx dt " f pxptqq verifies

dx ‹ dt " 0 " f px ‹ q (2.12)
stating that the state x ‹ is stationnary. In this case, defining a shifted system which dynamics is defined as fluctuations around x ‹ can be beneficial from a physical and numerical standpoint. This is the case for acoustical systems, exhibiting small fluctuations of pressure around the atmospheric pressure, or some electrical systems, like operational amplifiers that are supplied with a constant voltage source.

We now show how to establish a shifted finite dimensional pHs.

Consider a pHs of state variables x, Hamiltonian Hpxq and effort vector epxq " ∇ x Hpxq of form

(original pHs) ¨f w y ‹ ‹ ' F " S ¨e zpwq u ‹ ‹ ' E (2.13)
where the interconnection matrix S is constant.

We submit the system to a set of static inputs u " u ‹ under which we assume the existence of a state at rest x ‹ such that 9

x ‹ " 0. The corresponding pHs at rest reads

(pHs at rest) ¨0 w ‹ y ‹ ‹ ‹ ' F ‹ " S ¨e‹ zpw ‹ q u ‹ ‹ ‹ ' E ‹ (2.14)
with the constant outputs y ‹ , efforts e ‹ " epx ‹ q and memoryless variables w ‹ associated with e ‹ .

We define the shifted pHs by substracting Eq. (2.13) and (2.14) such as

(2.13) ´(2.14) ðñ F ´F‹ " S pE ´E‹ q (shifted pHs) ¨f w ỹ ‹ ‹ ' " S ¨ẽ zp wq ũ ‹ ‹ ' (2.15)
where we introduced the shifted: algebraic variables w " w ´w‹ ; memoryless laws zp wq " zp w `w‹ q ´zpw ‹ q; outputs ỹ " y ´y‹ ; input ũ " u ´u‹ and efforts ẽ " epxq ´e‹ . Finally, we define the shifted energy variables x " x ´x‹ , with

d x dt " dx dt ,

and the shifted Hamiltonian

Hp xq " Hp x `x‹ q ´e‹ ¨x ´Hpx ‹ q with e ‹ " ∇ x Hpx ‹ q, verifying ∇ x Hp xq " ẽ, (2.16)

where the second term of the right-hand side is tangent to the Hamiltonian at x " x ‹ (see dotted red line on Fig. 9) and will serve as the abscissa of the shifted Hamiltonian (see bottom plot). The last term is the energy of the system at rest, denoted by a red star. The parts of the Hamiltonian that are under the tangent (colored in yellow) will now be negative values of the shifted Hamiltonian Hpx ‹ q.

The shifting procedure summarises as follows.

Step 1: define the set of static inputs u ‹ ;

Step 2: compute the associated state at rest x ‹ that is solution of Eq. (2.14);

Step 3: define the shifted quantities in Eq. (2.15);

Step 4: compute the shifted state vector x, the shifted Hamiltonian and shifted efforts of Eq. (2.16);

Step 5: formulate the shifted pHs.

This procedure will be used on a model of the vocal tract (see chapter 3) to build signals that are fluctuations around a state at rest (as in acoustics) without neglecting nonlinearities.

Remark 4 (shifted pHs): Only Ref. [START_REF] Monshizadeh | Conditions on Shifted Passivity of Port-Hamiltonian Systems[END_REF] proposed a mathematical analysis of the method, tackling the conditions under which the system remains passive. To our knowledge, only our work of Refs. [WHS19; WHS20] applied this procedure to physical systems.

x = 0 x 0 H(x ) e(x ) H(x ) e(x )x H(x)
Hamiltonian H(x)

H(x) x = 0 0 Shifted Hamiltonian H(x) H(x)
Figure 9: Example of the shift procedure for a one-dimensional Hamiltonian function.

summary and perspectives

In this chapter, we introduced the port Hamiltonian framework, a powerful modelling tool that ensures the power balance through a conservative geometric structure encoding the power exchanges between the physical components. In the third section, we recalled and improved the discrete gradient (DG) numerical method to deal with cancelation errors that arise in certain situations. Finally, we presented a procedure to shift a given pHs around its state at rest to improve its physical interpretability and numerical conditioning.

The IDG method can be further improved. As of now, the switch between the finite difference and the alternative approximation of the discrete gradient operator could be smoothed out by using a linear interpolation between the value of the finite difference at δx k " DG and the alternate approximation. This could also improve the convergence of the Newton-Raphson implicit method.

The shifting procedure can be extended for the general case, where the interconnection matrix is not canonical, and generalised for the infinite dimensional case. This method can also be useful for the study of bifurcations that occur at a particular functioning point.

Part II

P O W E R -B A L A N C E D M O D E L L I N G O F T H E V O C A L T R A C T P O W E R -B A L A N C E D L U M P E D -PA R A M E T E R M O D E L L I N G O F T H E V O C A L T R A C T
The vocal tract is comprised of several cavities (see figure 10) surrounded by tissues.

To articulate vowels and consonants, the articulators (the muscles within the face and tongue) modify the geometry of the cavities (e.g. the oral cavity). Most models describe the vocal tract from an acoustical standpoint, either with high-dimensional models [VH Š08;[START_REF] Kunitoshi | Three-dimensional acoustic field in vocal-tract[END_REF] or with simpler approaches, like the lumpedparameter models [Mae82; ST95] or the digital ladder filters [KL62; VK94]. Recently, high-dimensional models allow for the use of monolithic FSI methods, e.g. in Refs. [START_REF] Guasch | A stabilized finite element method for the mixed wave equation in an ALE framework with application to diphthong production[END_REF][START_REF] Schickhofer | Compressible flow simulations of voiced speech using rigid vocal tract geometries acquired by MRI[END_REF], but this approach remains rare.

The modularity and scalability of lumped-parameter models have already allowed for many improvements on the first model of Flanagan [START_REF] Loton | Speech analysis synthesis and perception[END_REF] such as the addition of noise sources [START_REF] Birkholz | Construction And Control Of A Three-Dimensional Vocal Tract Model[END_REF], losses due to turbulences [START_REF] Birkholz | Simulation of losses due to turbulence in the time-varying vocal system[END_REF], a self-oscillating model of the tongue [START_REF] Elie | Extension of the single-matrix formulation of the vocal tract: Consideration of bilateral channels and connection of self-oscillating models of the vocal folds with a glottal chink[END_REF], and lateral channels in the mouth [START_REF] Elie | Simulating alveolar trills using a two-mass model of the tongue tip[END_REF] (for the synthesis of liquid consonants, e.g. [ ]).

None of the aforementioned models include guarantees regarding the energetic consistency of the models or simulations, or investigate the impact of the motion of the articulators on the flow. Usually, to account for the continuously changing geometry, driven by the articulators, while keeping an acoustical description, a common simplifying hypothesis is to neglect the energetic impact of the movements of the tissues (caused by the movement of the articulators) on the fluid, the so-called quasistatic assumption. Using this hypothesis leads to a paradoxical situation where the walls provide energy to the main flow, but the flow is modelled as if the walls were motionless (breaking the physical power-balance and passivity on the dynamics).

In this chapter, we propose a passive FSI model of the vocal tract that accounts for the transverse motion of the fluid induced by the displacement of the articulators in an energetically-consistent way. To do so, we approximate each of the time-varying waveguides (see solid colors on Fig. 10) by a FSI problem, as the one depicted on Fig. 11a. To build a simple and lightweight model of this system, we choose a lumpedparameter approach, in which each cavity is divided into subsystems called "tracts" (separated by dotted lines on figure 11a).

For each tract, we build a finite dimensional model using a set of geometrical approximations and physical hypotheses. Finally, we obtain the full vocal tract model by assembling tracts together (see Fig. 11b). Throughout this approach (modelling and assembly), the conservation of mass, momentum and energy is guaranteed by the use of the port-Hamiltonian (pH) framework. To connect the tracts, we consider basic junctions relying on the algebraic equations of mass balance and energy conservation (which together induce the momentum balance). At a fluid-fluid interface, we consider the mass flow q and its power conjugate quantity (the specific total enthalpy, denoted ψ and defined below) as stated on Fig. 12. Conversely, for a solid-fluid (or solid-solid) interaction, the connection variables are the resulting force exerted by the solid on the fluid (or the other solid) and the normal velocity of the interface.

S W T1 S1 v F T2 S2 v F T5 S5 v F q ψ v F Solid Fluid (b) Macroscopic component representation.
In this work, we first focus on the modelling of the behaviour of the fluid, then assemble the full vocal tract. Our modelling approach unfolds in two main steps.

step 1: from physics to the phs lumped-parameter model of a single tract -Assuming a set of geometrical and physical hypotheses, we build the macroscopic quantities describing the state of a tract. The corresponding dynamics and inputs-outputs are recast as a pHs, and discussed and interpreted using an equivalent electrical circuit representation. This first step is described in section 3.1

q 1 ψ 1 q 2 ψ 2 q K ψ K k ∈ {1, . . . , K} Sys 1 Sys 2 Sys K K k=1 q k = 0 ψ k = ψ k-1
Figure 12: Junction of K tracts interacting through a common interface (K " 2 for e.g. tracts T 1 and T 2 on Fig. 11a, K " 3 for e.g. the nasal derivation on Fig. 10). Note that at the interface, the power balance reads ř k q k ψ k " 0; the dual relation of the mass conservation is the continuity of the specific total enthalpy. step 2: assembling the full vocal tract -The connection of adjacent tracts with identical input-output configurations results in an incompatible connection setting. Therefore, tracts have to be interconnected through their energy (see p. 20 in chapter 2). In Sec. 3.2, we solve this problem in a way similar to Refs. [Mv92; Wu+14; Car16] by projecting the constrained pHs into a constrained subspace.

power-balanced lumped-parameter model of a single tract

In this section we propose a generic lumped-parameter model of a single tract, with a focus on the fluid flow. geometry of a tract -Each tract T is approximated by (H0) the spatial domain Ω T ptq " px, y, zq P r´ 0 , 0 s ˆr0, hptqs ˆr0, L 0 s ( depicted on Fig. 13. Its boundary BΩ T ptq is split into: the open and motionless surfaces S L and S R1 (at x "

Hypotheses and Approximations

v W e y F W e y ψ L q L e x ψ R q R e x -0 0 x = 0 y = 0 h(t) S R S L S W S B e y e x Ω T (t)
´ 0 and x " 0 , respectively); the impervious and uniformly moving top surface S W (located at y " hptq, moving with transverse velocity v W ); and the rigid motionless bottom surface S B (at y " 0). The volume of the tract Vptq is given by Vptq " 2 0 L 0 hptq.

Remark 5 (bottom surface): For the sake of simplicity, the bottom surface S B at y " 0 is assumed to be rigid and motionless. To adapt our work to the case where S B is a uniformly moving surface, use the methodology of Ref. [START_REF] Hélie | Self-oscillations of a vocal apparatus: A porthamiltonian formulation[END_REF].

physical hypotheses -We remind the reader of the physical hypotheses we consider (see chapter 1):

(H1) the fluid is inviscid (no viscosity); (H2) the flow is irrotational.

(H3) the fluid fullfils the perfect gas assumption;

(H4) the fluid undergoes an adiabatic and isentropic transformation (reversible process); (H5) the heat capacity at constant volume C V is assumed to be constant w.r.t.

the temperature within the considered range of fluctuations;

to which we add the lumped-parameter assumption (H6) and a linearising hypothesis (H7):

(H6) the volumetric mass density ρ is homogeneous within Ω;

(H7) the fluctuations of the mass density ρ are small compared to its value at the state of rest ρ 0 so that pρ´ρ 0 q {ρ 0 ąą 1.

The latter only impacts the thermodynamical behaviour of the fluid. (H7) enables us to consider a quasi-quadratic expression of the thermodynamical energy, better suited to a numerical context where the acoustical fluctuations of the mass density ρ ´ρ0 are small.

Within Ω T ptq, the local dynamics of the flow is governed by the momentum and mass conservation equations, i.e. :

Bρ Bt `divpρvq " 0, and Bρv Bt `divpρv b vq `gradppq " 0 (3.1)
where the volumetric mass density ρ, Eulerian velocity v and pressure2 p are physical fields, function of space and time.

boundary conditions -At the moving and impermeable boundary S W , the continuity of normal velocity relates the velocity of the fluid to the one of the boundary v B :

pv ¨nq ˇˇS W " v B ¨n which yields @x P r´ 0 , 0 s , v px, y " hptq, tq ¨ey " v W ptq, (3.2) where n " e y is the unit normal vector of the top surface. The power dual quantity of the (uniform) normal velocity of the wall v W is the force exerted by the wall onto the fluid F W . The product of both quantities is the power supplied to the flow through the moving interface

P ext,S W " F W v W .
The interactions at the open and motionless boundaries S j " tS L , S R u are described by the (extensive) mass flow q j and its conjugated power variable the (intensive) total specific enthalpy ψ j q j " ij S j ρv ¨ex dS and ψ j "

1 S j ij S j `1 2 v ¨v `p ρ dS @j P tL, Ru (3.3)
where denotes the specific internal energy, the expression of which is given in the next section. The mass flows are oriented according to the following convention: q L is the ingoing mass rate (in kg ¨s´1 ); q R is the outgoing mass rate (see the red arrows at the boundary on Fig. 13). The power flows at S L and S R read P ext,S L " ´qL ψ L and

P ext,S R " q R ψ R .
approximated velocity field and lumped mass density -Here, we supply approximations for the continuous velocity field v and the volumetric mass ρ. We generalise the kinematics proposed in ref [START_REF] Lopes | Energy Balanced Model of a Jet Interacting With a Brass Player's Lip[END_REF] to account for compressible phenomena so that, within each tract, we consider that vpx, y, tq « vpx, y, tq " ˜1 0 ¸vmx ptq Regarding the volumetric mass ρ, we use a (H6) piecewise constant approximation such as ρpx, y, tq « 1 Ω T ptq ¨ρptq, @px, yq P Ω T ptq.

(3.5)

Notation 1. For concision, the time dependence of v mx ,v y , v c , ρ and V and their derived quantities is omitted herafter.

Remark 6 (compression component of the velocity field): Computing the divergence of the volumetric momentum yields

divp ρ vq " ρv c 0 . (3.6)
When v c ą 0, the tract ejects matter, as the divergence is positive. Conversely, the tract pumps matter for negative values of v c . Thus, v c is directly related to the compression state of the fluid contained in the tract.

We now use these hypotheses and approximations to build the macroscopic quantities that describe the model, also called lumped-parameters.

Lumped-Parameter Model

energy of the flow -The total energy is the sum of the kinetic energy E kin and internal energy E comp . Integrating the kinetic energy density over the domain with the supplied fields v and ρ (eqs. (3.4) and eq. (3.5)) yields

E kin " ¡ Ω T ptq 1 2 ρv ¨v dV (3.4) « (3.5) 1 2 ρV « v 2 mx `1 3 v 2 y `1 3 ˆ 0 h v y ´vc ˙2ff . (3.7)
The internal energy E comp , hereafter called "compression energy", is obtained using hypotheses (H3-7) and field approximation (3.5) (see App. B.1 for the detailed computations) so that

E comp (H3-7) " P 0 V « γ 2 ˆρ ´ρ0 ρ 0 ˙2 ´1ff , (3.8)
where γ " Cp {C V and where P 0 is the pressure of the gas at its rest state (atmospheric pressure).

dynamics -The tract energy relies on 5 degrees of freedoms (DoFs): v mx , v c and v y for the kinematics; ρ and V for the "compression". The dynamics of the volume

-0 0 0 0 h(t) y v mx -0 0 0 x v y -0 0 0 v c

V is

d dt Vptq " S W d dt hptq " S W v W .
To compute the dynamics of the other DoF, we consider the weak form of the mass and momentum balances (3.1), i.e., projection onto a test function φpx, y, tq (see App. B.2). First, using φ " 1 Ω T ptq on the mass balance yields d dt ρ " ´vc 0 ρ.

(3.9)

We then test the momentum balance against functions 1 Ω T ptq e x and y e y so that dv mx dt " 1 2 0 pψ L ´ψR q and ρV 3

dv y dt " F W `SW p (3.10)
where p is the volume averaged pressure (see appendix B.3) and F W is the force applied by the wall unto the fluid. These equations simply state that the axial and transverse accelerations are a consequence of an enthalpy gradient, or a consequence of the effective wall force, respectively, without contributions related to the change of mass or geometry that would appear if momentums were considered as DoFs instead of the velocities (see Ref. [START_REF] Wetzel | Power balanced time-varying lumped parameter model of a vocal tract: modelling and simulation[END_REF]).

Finally we consider the test function φ " x e x . The resulting dynamics can be divided into 3 contributions: external compressive efforts (3.11a); the velocity of the wall modifying the volume, thus compressing the fluid (3.11b); the fluid resisting to the previous contributions (3.11c). The dynamics of v c is shown as such

V 3 dv c dt " ´S pψ L `ψR q ` 0 L 0 S F W (3.11a) `vy pS W vc {3 ´vyL0 2 0{hq (3.11b) ´2S " 1 2 `v2 mx `v2 c{3 `v2 y{3 ˘` p {ρ p 2 0{h 2 `1q ` p ρq  (3.11c)
where S " hL 0 .

causality of the boundary conditions -Looking at the equations of dynamics (3.10) and (3.11) yields the following causality. ψ L , ψ R and F W are the inputs of the system as they are unknown quantities, appearing as source terms in the dynamical equations (3.10) and (3.11). It follows that their power dual quantities q L , q R and v W " v y are the outputs which can be retrieved as q L " ρS rv mx ´vc `vy 0{hs , q R " ρS rv mx `vc ´vy 0{hs and v W " v y , (3.12) using the boundary conditions and proposed field approximations.

Casting the Lumped Model as a pHs

This subsection addresses the macroscopic port Hamiltonian formulation of a tract.

parametrisation of the hamiltonian -To make the interconnection matrix canonical 3 , we exhibit extensive energy variables pν L , ν R , Π y q ensuring trivial relationships between the internal efforts and the physical port variables pq L , q R , v W q of Fig. 11b. The state vector x T (size n x " 5) of a given tract is thus x T " " ν L ν R Π y m h ı with ν L :" 0 pv mx `vy 0{3h ´vc {3q , ν R :" 0 pv mx ´vy 0{3h `vc {3q (3.13a)

Π y :" v y ρV 3 , m " ρV and h " V S W , (3.13b)
which can be split into a set of axial kinematic variables ν L and ν R , a transverse incompressible component Π y , and a set of compression variables x c " rm, hs . The Hamiltonian H T px T q is the sum of the kinetic and compression energies expressed as functions of the energy variables x T such as:

H T px T " " x kin ν L ν R Π y m h xc ı q " H kin px kin q `Hcomp px c q " m 2 2 0 `ν2 L `ν2 R ´νL ν R ˘`3Π 2 y 2m `P0 hS W « γ 2 ˆρpm, hq ´ρ0 ρ 0 ˙2 ´1ff (3.14)
where ρpm, hq " m {pS W hq " ρ is a function that reconstructs the volumetric density.

Notation 2. The function ρ enables us to keep the energy in an acoustical friendly form, where ρ would be a natural choice of energy variable. In the following and for clarity, we only specify the dependence of ρ on m and h when needed.

With the given choice of energy variables, the effort vector e T reads

e T " » - - - - - - - - e ν L " q L e ν R " q R e Πy " v y e m " ψ Ω e h " ´ p Ω S W fi ffi ffi ffi ffi ffi ffi ffi fl " » - - - - - - - - - p2ν L ´νR q m {2 2 0 p2ν R ´νL q m {2 2 0 3Πy {m pν 2 L `ν2 R ´νL ν Rq 2 2 0 ´3Π 2 y 2m 2 `P0 γ ρ 2 0 pρpm, hq ´ρ0 q S W h P 0 S W `γ 2 ´1˘´P 0 γ 2 ´ρ ρ 0 ¯2 S W fi ffi ffi ffi ffi ffi ffi ffi ffi fl (3.15)
where each effort has a physical meaning: by design, the first three are the mass rates q L and q R and the wall velocity; e m is the volume-averaged specific total enthalpy ψ Ω ; and e h is the reaction to the net pressure force of the fluid on the wall (see Eq. (B.1.5)).

dynamics -Using Eq. (B.2.9) with φ " 1, the dynamics of the mass yields 9 m " q L ´qR " e ν L ´eν R stating that the mass rate 9 m equals the sum of the incoming mass flows.

The dynamics of the height reads 9 h " v W " v y " e Πy . The dynamics of Π y is derived by using Eq. (B.2.10) with φ " y {h e y which yields

9 Π y " 9 m Π y m `FW ` p S W " pe ν L ´eν R q Π y m ´eh `FW . (3.16)
The dynamics of ν L and ν R is obtained by using Eq. (B.2.11) with the test functions φ " 0{m p1 ´x{ 0 q e x for ν L and φ " 0{m p1 `x{ 0 q e x for ν R4 which leads to

d dt ν L " ψ L ´Πy m v y ´em and d dt ν R " ´ψR `Πy m v y `em .
(3.17) interconnection matrix -We identified the components of the effort vector e T in the equations of dynamics. However, some of them contain a term modulated by Πy {m, thus requiring the introduction of a gyrator [Lop16, Def. 8, p. 24] to make the interconnection matrix canonical. This gyrator is modelled as a conservative memoryless component converting the mass accumulation 9 m and the transverse velocity v y into the enthalpy vy Πy {m and the force Πym {m. Let us consider the vectors w and zpwq (of size n w " 2) so that

w T " « 9 m " e ν L ´eν R v y " e Πy ff , zpw T q " Π y m « 0 1 ´1 0 ff ¨wT " ˜vy Πy {m ´9 m Πy {m ¸, (3.18)
verifying zpw T q w T " 0.

The input vector u T gathers the inputs ψ L ,ψ R and F W so that u T " rψ L , ψ R , F W s . Similarily, the output vector is defined as y T " rq L , ´qR , v W s .

Then, the full interconnection matrix of the pHs reads

¨fT w T ´yT ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' " ¨9 ν L 9 ν R 9 Π y 9 m 9 h 9 m v W ´qL q R ´vW ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' " » - - - - - - - - - - - - - - - - - - - - - - ´1 ´1 `1 `1 `1 ´1 ´1 ´1 `1 `1 ´1 `1 `1 ´1 `1 ´1 `1 ´1 fi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi fl S T ¨eν L " q L e ν R " q R e Πy " v y e m " ψ Ω e h " ´ p Ω S W v y Πy m ´Πy m 9 m ψ L ψ R F W ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' (3.19
) where the matrices J, G w and G p can be identified according to the structure presented in Eq. (2.4). As the full interconnection matrix is canonical, S T will be identical for each tract in the vocal tract network.

Electrical Equivalent Circuits

While bond graph would be a possible graphical representation for the interactions of the components, we here choose the electrical equivalent circuits for their sim- 1) between mechanical and electrical quantities [START_REF] Floyd A Firestone | A new analogy between mechanical and electrical systems[END_REF] and laws from which one can build an equivalent electrical circuit. This graphical representation highlights the interactions between the different components, thus allowing for physical interpretation. From a pHs standpoint, the interconnection matrice J acts as the PCB of the circuit as it only encodes the connections between the components.

In our model, the axial and transverse flows involve different boundary flows and efforts variables. The simplest electrical equivalent representation is obtained by using two different analogies (see table 1): the fluid mechanical one for the axial flow; and the solid mechanical one for the transverse motion of the fluid flow and of the wall.

The resulting equivalent electrical circuit of the tract model is depicted on Fig. 15. Note that all elements (gyrator, inductors and capacitors) have a non linear behaviour as their attributes (conversion factor, inductance and capacitance) vary with the state of the system. (axial flow to transverse motion). It is also responsible for the change of analogy between the axial and transverse flows. The subsystem comprised of the gyrator and both capacitors acts as a gas spring accounting for the compressibility of the fluid in a time varying volume.

The following degenerated study cases provide further insight on the coupling role of the gyrator: Case 1 (motionless top surface, v y " v W " 0): the top surface is assumed to be motionless with v W " v y " 0 so that the voltage drop of the bottom branch of the gyrator vy Πy m vanishes, connecting the capacitor m to the ground. Compression then corresponds to mass change only, without any transverse motion.

Case 2 (steady mass flow, q L " q R ): this implies a null total mass rate 9 m " 0. The voltage across the top branch of the gyrator now reads 9 mΠy m " 0 which directly connects the transverse capacitor to the ground. Compression of the fluid thus translates into volume change only, with a constant fluid mass in the tract. The most simple application of this case is to consider q L " q R " 0, meaning that the tract acts as a piston.

We briefly compare our model the popular one of Maeda [START_REF] Maeda | A digital simulation method of the vocal-tract system[END_REF], the equivalent electrical circuit of which is shown on Fig. 16. In both cases, the circuit has the characteristic T-cell design, with inductors on the series branch and the capacitor on the parallel branch. The model of Maeda is linear and does not use any coupling element at the junction. Therefore, the mass flow induced by the movement of the wall is in parallel with the compression effect.

u L R f L f u L L f u R R f u R L w u w R w C w u s C f u c
The upper branch of our model (inductor h and inductor Π y ) is not present on the model of Maeda, meaning that the model does not account for the transverse motion and volume change. Finally, Maeda added a basic fluid dissipation model5 to enhance the physical realism, especially in the case where the cross section is small.

Shifted Tract Model

The model we established is a fluid mechanical model. To achieve a (nonlinear) acoustical representation, the physical variables have to be expressed as fluctuations around a state at rest that usually corresponds to atmospherical conditions. This procedure also leads to improvements of the numerical performances.

To do so, we use the shift procedure introduced in chapter 2, within Sec. 2.3.4. We assume that the system (3.19) is under atmospheric conditions, meaning that the surrounding air is at rest at atmospheric pressure P 0 . Thus, the inputs reads

u ‹ " " ψ ‹ L ψ ‹ R F ‹ W ı " " 0 0 ´P0 S W ı under
" ν ‹ L ν ‹ R Π ‹ y m ‹ h ‹ ı " " 0 0 0 m 0 h 0 ı
. An alternative and more complex choice would be to shift the system in the vicinity of the trajectory of hptq and mptq but is not addressed here. At the rest state, the effort vector reads e ‹ T " r0, 0, 0, 0, ´P0 S W s and the output vector is y ‹ " r0, 0 0s . The latter corresponds to a uniform atmospheric pressure in the tract exhibiting no power exchanges with the exterior.

The shifted energy variables are gathered in xT " " νL νR Πy m hı and the shifted Hamiltonian Hp xT q reads Hp xT q " Hp xT

`x‹ T q ´eT px ‹ T q ¨x T ´Hpx ‹ T q " p m `m0 q 2 2 0 `ν 2 L `ν 2 R ´ν L νR ˘`3 Π2 y 2p m `m0 q `P0 S W γ 2 ˆρp m, hq ρ 0 ˙2 p h `h0 q
where ρp m, hq " ρp m `m‹ , h `h‹ q ´ρ0 is a function that reconstructs the fluctuations of the volumetric mass around its state of rest value ρ 0 .

The (shifted) efforts ẽT are ẽT "

¨p2 νL ´ν R q p m`m 0 q {2 2 0 p2 νR ´ν L q p m`m 0 q {2 2 0 3 Πy {p m`m 0 q p ν2 L `ν 2 R ´ν L νRq 2 2 0 ´3Π 2 y 2p m`m 0 q 2 `P0 γ ρ 2 0 ρp m, hq ´P0 γ 2 ´ρ ρ 0 ¯2 S W ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ '
Now the coordinates xT " 0 coincide with the state at rest of the system (under atmospheric conditions). The interconnexion matrix is not altered by the shift procedure as it canonical (no dependence on the physical parameters and energy variables6 ). The algebraic-differential formulation then reads

¨B xT Bt wT ´ỹ T ‹ ‹ ‹ ' " S T ¨ẽ T zT p wT q ũT ‹ ‹ ‹ '
, where we introduced the (shifted) inputs ũT and outputs ỹT ũT "

" ψL " ψ L , ψR " ψ R , FW " F W ´P0 S W ı , ỹT " " qL " q L , qR " q R , ṽw " v W ı ,
and shifted variables wT and algebraic laws zT p wT q wT " « 9 m " qL ´q R ṽy ff and zpw T q " Πy m `m‹

« 0 1 ´1 0 ff wT .

Summary

In this section, we derived a two-dimensional power-balanced lumped-parameter model of a single tract. It accounts for the transverse flow induced by the displacement of the articulators. Using a set of physical and geometrical hypotheses and approximations, we discretised the associated PDEs and energy function to produce a macroscopic model. We then cast it as a canonical pHs, which allowed for its physical interpretation by the means of an equivalent electrical circuit. We highlighted the presence of a new coupling element, a gyrator, that is responsible for the reciprocal power transfers between the axial and transverse flows. We briefly compared our model to the well-know model of Maeda [START_REF] Maeda | A digital simulation method of the vocal-tract system[END_REF]. Finally, we established an equivalent shifted pH formulation of the tract model to achieve an acoustical description.

To build the vocal tract, we now assemble multiple tract models together.

assembly of the full vocal tract

The section addresses the assembly of several tracts T k Pt1,...,Ntu through conservative junctions η jPt1,...,N λ u as shown on Fig. 17. As stated in the introduction, we assume that at a junction η j

η1 η2 η3 η4 η5 η6 η7 η8 T 1 T 2 T 3 T 4 T 5 T 6 T 7 T 8 T 9 T 10 T 11 T 12
N ÿ k"1 d j,k q k " 0 with d j,k " $ ' ' ' & ' ' ' % `1 if q k is entering η j ´1 if q k is leaving η j 0 if q k is disconnected from η j (directed mass flow balance) (3.1a) and ψ 1 " ψ 2 " . . . " ψ N (enthalpy continuity) (3.1b)
involving (at the junction) N directed mass flows d j,k q k and enthalpies ψ k , @k P t1, . . . , Nu. The coefficients d j,k are extracted from the incidence matrix of the graph composed of the nodes η j and branches T k equipped with their mass flow orientation.

Note that connecting tracts together in this manner is strictly equivalent to connecting the electrical equivalent circuits (see Fig. 15) by a node, as Eq. (3.1) is the fluid mechanical analog of Kirchoff's current law and continuity of potentials (see Tab. 1). For instance, consider the junction η 1 on Fig. 17 where N " 4 so the mass balance reads q R1 " q L2 `qL3 `qL4 .

Considering more complex junction models would introduce additional storage and/or memoryless components at the interconnection. We propose such a component in chapter 5.

Problem Statement and Proposed Solution

At both its open boundaries, a tract is endowed with an input-output configuration imposed by the physics at hand. At a boundary, the input-ouput configuration is said to be enthalpy driven if the enthalpy ψ is the input, while it is said to be mass flow driven if the mass flow q is the input. The compatible connection of two tracts is usually ensured by the connection of complementary configurations (e.g. enthalpy driven configuration connected to a mass flow driven configuration). However, due to its symetrical form, our tract model only allows for enthalpy driven configurations at both extremities. This creates an incompatible connection (see chapter 2), as two contiguous tracts exhibit identical input-output configurations. This is also called a causality issue, which manifests itself as a constraint.

As a simple example, consider the two-tract example depicted on Fig. 18 where T 1 interacts with T 2 through its right interface S R1 (resp. S L2 , see equations below).

q R 1 = q L 2 T 1 T 2 {h 1 , m 1 , Π y1 } {h 2 , m 2 , Π y2 } v W 1 e y F W 1 e y v W 2 e y F W 2 e y ψ L1 q L1 e x ψ R2 q R2 e x Figure 18: A basic two-tract example. » - - - - - - - - - f T 1 w T 1 ´qL1 q R1 v W 1 fi ffi ffi ffi ffi ffi ffi ffi ffi fl " S T 1 ¨» - - - - - - - - e T 1 zpw T 1 q ψ L1 ψ R1 F W fi ffi ffi ffi ffi ffi ffi ffi fl Eq. : pHs of tract T 1 » - - - - - - - - - f T 2 w T 2 ´qL2 q R2 v W 2 fi ffi ffi ffi ffi ffi ffi ffi ffi fl " S T 2 » - - - - - - - - e T 2 zpw T 2 q ψ L2 ψ R2 F W fi ffi ffi ffi ffi ffi ffi ffi fl
Eq. : pHs of tract T 2 At the shared interface, the mass conservation reads q R1 ´qL2 " 0 where each mass flow results from the internal dynamic of each tract, thus making this rela-tion implicit. To connect the two tracts, we introduce the Lagrange multiplier ψ λ , constraining the dynamics of T 1 and T 2 such as Eqs. (3.1) hold:

¨0 ψ R1 ´ψL2 ‹ ‹ ' " » - - - `1 `1 ´1 ´1 fi ffi ffi fl ¨ψλ `qR1 ´qL2 ‹ ‹ ‹ ' . (3.2)
Remark 7 (electrical equivalent of the causality issue): Looking at the connection of two equivalent electrical circuits of Fig. 15, this correponds to the case where two inductors are in series and need to be merged into an equivalent inductors.

In Fig. 17, each junction η j correponds to a similar constraint. We now show how to properly formulate and solve this constraint problem.

To include the constraints in a pH formulation, we use the differential-algebraic port-Hamiltonian system (DA-pHs) formulation, introduced in chapter 2, the algebraicdifferential formulation of which we recall below

¨f w 0 y ‹ ‹ ‹ ‹ ' " » - - - - - - - J G w C λ G p ´G w J ww 0 G wp C λ 0 0 0 ´G p ´G wp 0 J pp fi ffi ffi ffi ffi ffi ffi fl ¨e zpwq λ u ‹ ‹ ‹ ‹ ' . (3.3)
A DA-pHs accounts for a set of N λ constraints with N λ Lagrange multipliers λ and a set of N λ algebraic equations relating the internal efforts e through a constraint matrix C λ (size n x ˆNλ ) such as 0 " C λ e .

For the two-tract example, the matrix

C λ " " e ν L1 e ν R1 e Π y1 e m 1 e h 1 e ν L2 e ν R2 e Π y2 e m 2 e h 2 0 1 0 0 0 ´1 0 0 0 0
ı selects the relevant quantities (q R1 and q L2 ) from the collected effort variables e " " e T 1 , e T 2 ı , thus encoding the mass balance. The vector of Lagrange multipliers reads λ " rψ λ s (see Eq. (3.2)).

As constrained systems might be hard to interpret, analyse and simulate, they are often recast into a regular ODE. One popular method is the index reduction by differentiation [START_REF] Kunkel | Differential-algebraic equations: analysis and numerical solution[END_REF], where the constraints equations are derived with respect to time to transform the algebraic relations into a set of regular ODE. However, it does not ensure that the system remains physical as "discretisations and round off errors may lead to numerical results that violate the constraints" [KM06, p. 273]. Moreover, replacing an algebraic equation by an ODE can be harmful as the problem can become stiff and require a very robust numerical method.

In a very similar modelling setting, Mora et al. [START_REF] Alejandro | Port-Hamiltonian Modelling of Fluid-structure Interactions in a Longitudinal Domain[END_REF] introduced a purely compressible cavity between tracts that acts as a causality-compatible junction and solves the causality issue beforehand. However, the physical interpretation of this junction is hard to make, while its role in numerical simulations is not negligible.

Other methods try to solve constraints by coupling the constrained components through their energy. We choose this class of methods because they are physically interpretable and ensure the passivity of the final model.

The equivalent component method [START_REF] Najnudel | Analog circuits and Port-Hamiltonian realizability issues: a resolution method for simulations via equivalent components[END_REF] is one of them. In our case, it corresponds to computing the equivalent inductor at a junction, which will be parametrised by the physical parameters of the contiguous tracts. We applied this method in Ref. [START_REF] Wetzel | Power balanced time-varying lumped parameter model of a vocal tract: modelling and simulation[END_REF], which enables us to connect two tracts together. However, this approach is not adapted for the connection of more than three tracts, as the nonlinear nature of the inductors and the propagation of the parameters result in a very complex combinatorics.

Thus, we consider the coordinate change method [Car16; Mv92; Wu+14; van13], a lesser-known method 7 , that casts a DA-pHs

(3.3) into a semi-explicit DAE of index 8 1 of form $ ' & ' % dχ θ dt " f pχ θ , χ λ q 0 " gpχ θ , χ λ q (3.4)
where the constrained states χ λ are determined by the implicit relation on the instanteneous (algebraic) part (hence the "semi-explicit" name), and act as parameters for the dynamics of the unconstrained states χ θ .

To address the causality problems of a complete network of tracts, we propose a systematic procedure that: (i) automatically casts a network of systems like the one of Fig. 17

From Directed Graphs to DA-pHs

This procedure is inspired by Ref. [START_REF] Müller | Fully-implicit Algebro-differential Parametrization of Circuits[END_REF]. For a more general context, see Ref. ??.

To formulate the DA-pH of a network of pHs, we use a directed graph composed of nodes η j as vertex and systems k as edges endowed with a flow orientation. Under analogies shown in Tab. 1, the proposed procedure is valid for a class of pHs detailed hereafter.

Consider N s pHs of form (2.3) (indexed by the variable k P t1, . . . , N s u " 1, N s ), of energy variable x k , and rearrange their inputs and outputs as

u k " « u λk ûk ff and y k " « y λk ŷk ff " « G λk Ĝ pk ff e k @k P 1, N s (3.5)
7 Which can be seen as a generalisation of the equivalent component method. 8 See [START_REF] Volker Mehrmann | Index Concepts for Differential-algebraic Equations[END_REF] for an overview of the index concept in the DAE theory where y λ k (resp. u λk ) is a vector gathering two constrained (coupled) outputs (resp. inputs) and where p ŷ, ûq denotes the remaining (unconstrained) external ports9 . Both G λk and Ĝ pk are obtained by rearranging the lines of G pk .

Graph-based Assembly

Consider a directed acyclic graph GpN , Eq (directed graph without self-loops [START_REF] Newman | Networks[END_REF]) with nodes N " tη 1 , . . . , η N λ u and edges E " t 1 , . . . , Ns u. As in Ref.

[Le 15], the j-th vertex represents the junction η j and the k-th edge k represents the k-th system.

Under the definition of Eq. (3.5), the tail of the k-th edge corresponds to the first component of the constrained input u λk (resp. output y λk ) and the head of the edge to the second component (see Fig 19).

η 0 η 1 η 2 η 3 η 4 η 5 η 6 η 7 η 8 η 9 η 10 η 11 η 12 T 1 T 2 T 3 T 4 T 5 T 6 T 7 T 8 T 9 T 10 T 11 T 12
Figure 19: Oriented graph of the network in Fig. 17; N s " 12 and N λ " 13.

As in Eq. (3.1a), the incidence matrix DpGq of G encodes the relations between the directed edges (systems) and vertex (nodes/junctions) (3.7)

rDpGqs jk " d jk " $ ' ' ' & ' ' ' % `1 if k is incident to η j ´1 if k leaves node η j 0 if k is disconnected from η j @k P 1,
where its i-th element is denoted by Cpiq. Note that N C can be empty.

Then, given a directed graph G, a set of pHs like in Eq. (3.5) and a set of control nodes N C , the DA-pH formulation of the constrained network is obtained algorithmically, as shown hereafter.

From directed graph to the DA-pH formulation

We generate the DA-pHs formulation of the constrained system as follows.

Build the concatenated energy variable vector X " The equation of dynamics can be formulated as

» - - - - - F W 0 ´Yfi ffi ffi ffi ffi fl " » - - - - - - - J G w C λ pGq G p ´G w 0 0 0 C λ pGq 0 0 0 ´G p 0 0 0 fi ffi ffi ffi ffi ffi ffi fl ¨» - - - - - E ZpWq λ U fi ffi ffi ffi ffi fl (3.8)
where F " 9 X, J " diag pJ 1 , . . . , J Nt q, G w " diag pG w1 , . . . , G wNt q and the external ports vectors U, Y and the matrix G p are where diag pA 1 , . . . , A p q is the diagonal block-matrix with blocks A 1 , . . . , A p .

U " » - - - - - - u C û1 . . .
In this DA-pH formulation, some objects are left to be defined: C λ pGq in step 1, λ in step 2, y C and u C in step 3 and ĜC pGq in step 4.

step 1: We now define the constraint matrix. Under Eq. (3.5), the constraint equation (3.1a) at a junction η j , j P N takes the following form

0 " Q j in `Qj out 0 " ÿ k"t1...Nsu d j,k "1 " 0 1 ı y λk `ÿ k"t1...Nsu d j,k "´1 " 1 0 ı y λk 0 (3.5) " ÿ k"t1...Nsu d j,k "1 " 0 1 ı G λk e k px k q `ÿ k"t1...Nsu d j,k "´1 " 1 0 ı G λk e k px k q 0 " C λ pGq EpXq (3.10)
where the selection vectors " 1 0 ı and " 0 1 ı extracts the constrained output associated with Q in or Q out (sum of incoming/outgoing flow). Note that if the k-th system is not connected to η j (d j,k " 0), then the selection vector is implicitly " 0 0 ı . In the case of the pH tract model given in Sec. 3.1.3, y λk " " q Lk , q Rk ı , " 0 1 ı y λk corresponds to ´qRk and " 1 0

ı y λk to q Lk .
Then, according to Eqs. (3.10) and (3.8), the constraint matrix C λ pGq (size pN λ ŃC q ˆpN s ř k n xk q) is defined by assembling the matrices G λk as

rC λ pGq s jk " c jk " $ ' ' ' & ' ' ' % r0 , 1s G λk if d jk " 1, r1 , 0s G λk if d jk " ´1, 0 otherwise, @j P I z C @k P 1, N s
where only the subset of nodes N z N C is spanned (not accounting for control nodes).

The number of constraints is now N λ " N λ ´NC .

step 2: Associate a Lagrange multiplier λ j to each constraint and gather them in the Lagrange multiplier vector λ such as " λ ı j " λ j @j P I z C.

step 3: According to Eq. (3.9), span the set of control nodes N C to define the vector of N C control outputs y C as

ry C s " ÿ k"t1...Nsu d Cp q,k "1 " 0 1 ı y λk `ÿ k"t1...Nsu d Cp q,k "´1 " 1 0 ı y λk
The vector of N C control inputs u C is defined in the same manner, substituting

y λk by u λk . If N C " H, then y C " u C " H.
step 4: The control port interconnection matrix ĜC (see Eq. (3.9)) is a N C ˆNs block matrix defined similarily as C λ , but on N C . Its entries g k , of size 1 ˆnxk , are

" Ĝ C ‰ k " g k " $ ' ' ' & ' ' ' % r0 , 1s G λk if d Cp qk " 1 r1 , 0s G λk if d Cp qk " ´1 r0 , 0s G λk if d Cp qk " 0 @ P 1, N C k P 1, N s
where we are only looking at the nodes that belong to the subset N C and the edges that are connected to these ( d Cp qk ‰ 0). If N C " H, then ĜC is empty (zero rows).

From DA-pHs to Semi-explicit pHs

We showed how to automatically write the DA-pH formulation of the network of systems that are coupled by constraints. We now recall the coordinate change used to derive the semi-explicit formulation of the above implicit DA-pH. As stated in Sec. 3.2.2.1, we now cast the DA-pHs into a semi-explicit form like in Eq. (3.15) below, using the method proposed in Ref. [START_REF] Van Der Schaft | On the Hamiltonian formulation of nonholonomic mechanical systems[END_REF][START_REF] Van Der Schaft | Port-Hamiltonian Differential-Algebraic Systems[END_REF].

step 1: constructing the coordinate change -Denote C K λ the left annihilator of C λ a matrix of size pn x ´Nλ q ˆnx that verifies C K λ C λ " 0. The coordinate change matrix P (size n x ˆnx ) is given by

P :" « C K λ pC λ C λ q ´1C λ ff (3.11)
where the second block of the matrix is the left Moore-Penrose inverse of C λ , as proposed in Ref. [START_REF] Wu | Port Hamiltonian System in Descriptor Form for Balanced Reduction: Application to a Nanotweezer[END_REF]. Note that the use of the Moore-Penrose inverse is only valid if C λ does not depend on the physical parameters of the system.

step 2: new parametrisation of the hamiltonian -The new state vector χ is divided accordingly to the structure of P such as

χ " Px " » - χ θ χ λ fi fl and f 1 " dχ dt " « f θ f λ ff (3.12)
where we introduced the unconstrained energy variables χ θ , of size n θ " n x ´Nλ , the constrained energy variables χ λ , of size N λ , and their corresponding fluxes, f θ and f λ . This dichotomy will be applied to other objects throughout the method.

Let us define the corresponding Hamiltonian H 1 pχq and efforts variables e 1 as H 1 pχq " Hpx " P ´1χq and e 1 " P ´ e " ∇ χ H 1 pχq " . The term PC λ λ thus only contributes to f λ , which means that the dynamics 9 χ θ does not depend on λ anymore. Moreover, the constraints read ´pPC λ q e 1 " 0 trims down to e λ pχ θ , χ λ q " 0. This last equality enables the implicit evaluation of χ λ as a function of χ θ (instead of integrating the dynamics 9 χ λ ).

« e θ
As a consequence, the Lagrange multipliers λ can be removed from the DA-pHs 10 (but may still be recovered from the dynamics 9 χ λ ). The splitting and removal of 9 χ λ and λ lead to the following semi-explicit form: with J θ pχq, G θ w and G θ p defined after the splitting of the matrices

PJpxqP " « J θ J θλ ´J θλ J λ ff , PG w " « G θ w G λ w ff and PG p " « G θ p G λ p ff . (3.16)
In this section, we have shown how to formulate the full network of pHs as a semi-explicit system using directed graphs and a coordinate change. We now apply this method to the case of the vocal tract.

Application to the Vocal Tract

A blind application of the above assembly method is not optimal, especially when dealing with large networks of systems. Therefore, we vary the assembly procedure by reordering variables using a permutation matrix to simplify the computations and the final system. This enables us to exploit the particular structure of the tract model.

Notation 3. For legibility, we use the notation for the original tract model, before the shift (see Sec. 3.1.3). As the interconnection matrix is canonical, the assembly method applies identically in both cases (shifted and non-shifted).

DA-pH formulation of the vocal tract

Consider a set of N t tracts, the assembly configuration of which is given by the graph G T . The full state vector of the vocal tract X VT is defined by setting 11 x k " x T k (size n xk " 5) and using the permutation matrix B x (size N t ¨nxk ˆNt ¨nxk ) which gives where X VT is split into: the axial flow inertial components x ax , the transverse motion inertial components x tr and the compression components x comp " rx m , x h s . The Hamiltonian of the full vocal tract is

X VT " B x » - - - x T 1 . . . x T N t fi ffi ffi fl " » - - - - - - - " ν L1 , ν R1 , . . . ν LNt , ν RNt ı " Π 1 , . . . , Π Nt ı " m 1 , . . . , m Nt ı " h 1 , . . . ,
H V T pX VT " "
x ax x tr x m x h ı q " H ax px ax , x m q `Htr px tr , x m q `Hcomp px comp q sorted by axial flow, transverse motion and compression energies.

To facilitate the assembly process, we take advantage of the quadratic nature of the kinetic Hamiltonian w.r.t. the kinetic energy variables ν L , ν R and Π y .

For the full axial flow, we get H ax px ax , x m q " 1 2 x ax Q ax px m qx ax with

Q ax px m q " diag pµ 1 , . . . , µ Nt q and µ k " m k 2 0k « 1 ´1 2 ´1 2 1 ff . (3.18)
The total energy of the transverse motions is

H tr px tr , x m q " 1 2 x tr Q tr px m qx tr with Q tr px m q " 3 diag px m q ´1.
The total compression energy reads

H c ´xcomp " " m 1 . . . m Nt h 1 . . . , h Nt ı ¯" Nt ÿ k"1 H comp px ck " rm k , h k sq .
Accounting for the permutations, the efforts read

E VT " B x » - - -
« Q ax px m qx ax Q tr px m qx tr ff " » -´qL1 q R1 q L2 . . . q LNt q RNt ¯ ´vy1 v y2 . . . v yNt ¯ fi fl (3.19)
thanks to the permutation matrix B x .

We direct each tract T k so that the direction of the edge coincides with the main direction of the flow (from left to right, see Fig. 13). Therefore, the head (resp. the tail) of edge k correponds to the right surface S Rk (resp. the left surface S Lk ) so that G λk and Ĝpk are

G λk " « 1 0 0 0 0 0 ´1 0 0 0 ff and Ĝpk " " 0 0 1 0 0 ı .
For the memoryless components, we choose w k " w T k and zpw k q " zpw T k q (see Eq. (3.1.5)). The full DA-pH formulation of the constrained vocal tract is

» - - - - - F VT W VT 0 ´Y fi ffi ffi ffi ffi fl " » - - - - - - - J VT G w,VT C λ pG T q G p,VT
´G w,VT 0 0 0

C λ pG T q 0 0 0 ´G p,VT 0 0 0 fi ffi ffi ffi ffi ffi ffi fl ¨» - - - - - E VT Z VT pW VT q λ U fi ffi ffi ffi ffi fl (3.20)
where the interconnexion matrix J VT is obtained by applying the permutation matrix B x to the concatenated interconnexion matrices of each individual tract J k such as

J VT " B x diag pJ 1 , . . . , J Nt qB x .
Similarily, for the gyrators we have

G w,VT " diag pG w1 , . . . , G wNt qB ´1 x .
example 1: two coupled tracts -The directed graph that corresponds to the two tract model of Fig. 18 is shown on Fig. 20 where

N C " tη 1 , η 3 u implying N λ " 1. η 1 η 2 η 3 T 1 T 2 q R1 = q L2
Figure 20: Directed graph of the two-tract example (equivalent to Fig. 18).

Therefore, the constraint matrix gives C λ " " 0 1 ´1 0 0 0 0 0 0 0 ı , and G p,VT reads

G p,VT " » - - - - - - - q L1 q R1 q L2 q R2 v W 1 v W 2 e m 1 e m 2 e h 1 e h 2 `1 0 0 0 0 0 0 ´1 `1 0 0 1 O p2ˆ2q O p2ˆ4q O p2ˆ4q O p2ˆ4q fi ffi ffi ffi ffi ffi ffi fl encoding Y VT " G p,VT E VT " " ´qL1 q R2 v W 1 v W 2 ı
where the external ports q L1 and q R2 are control ports for the main axial flow (instead of being constrained). The horizontal lines in the vectors/matrices are related to the notations of Sec. 3.2.2. Above the matrix, we indicate to which component of the effort vector the columns of G p,VT are linked.

q R1 = q L2 q R2 = q L3 q R3 = q L4 q R4 = q L5 ν 12 ν 23 ν 34 ν 45 T 1 T 2 T 3 T 4 T 5 {h1, m1, Πy1} {h2, m2, Πy2} {h3, m3, Πy3} {h4, m4, Πy4} {h5, m5, Πy5} vW 1 FW 1 vW 2 FW 2 vW 3 FW 3 vW 4 FW 4 vW 5 FW 5 ψ L1 -q L1 ψ R5 q R5
Figure 21: Example of a vocal tract with N t " 5. In green, the assembly constraints at each junction. For each tract, the energy variables not affected by the assembly are m k , h k , Π yk @k P t1, . . . , N t u. The equivalent energy variable at each junction is

ν k,k`1 @1 ě k ě N t ´1
example 2: tube with variable cross section -We now look at a tube composed of N t tracts in series (depicted on Fig. 21 for N t " 5 tracts), the directed graph of which is shown on Fig. 22, containing N λ " N t ´1 constraints. Setting

η 1 η 2 . . . η k . . . η Nt ηN t+1 T 1 T 2 T k-1 T k T Nt-1 T Nt q R1 = q L2 q Rk-1 = q Lk q RN t-1 = q LN t
N C " tη 1 , η Nt`1 u (with N C " 2) gives that N λ " N t ´1 ´NC " N t ´3 such as C λ (of size N λ ˆnx ) reads C λ " » - - - - - q L1 q R1 q L2 q R2 q L3 q RpNt´1q q LNt q RNt e Π e m e h `1 ´1 `1 ´1 `1 ´1 O pN λ ´NC ˆ3Ntq fi ffi ffi ffi ffi fl encoding q R1 ´qL2 " 0 q R2 ´qL3 " 0 . . . q RpNt´1q ´qLNt " 0 .
The matrices connecting the efforts to the output ports reads

G p,VT " » - - - - - q L1 q R1 q L2 q RNt 1 0 0 ´1 e Π e m e h O p2ˆp2Nt´2qq O p2ˆNtq O p2ˆ2Ntq O pNtˆ2Ntq I pNtˆNtq O pNtˆ2Ntq fi ffi ffi ffi ffi fl so that Y VT " ´G p,VT E VT " " q L1 , ´qRNt v W 1 , . . . , v W Nt ı .

Semi-explicit Formulation of the Vocal Tract

We now apply the coordinate change on the DA-pH of the vocal tract. Utilizing the permutations and the previous matricial definitions of the Hamiltonian yields a simplified final model, thus indicating that blindly applying the method is not the most effective way.

Here, only the axial kinematic efforts variables are constrained. Accounting for this and the permutations, P simplifies as

P " » - - - - - - - - - - - - - x ax x tr x m x h P ν O p2Nt´N λ ˆ3Ntq O p3Ntˆ2Ntq I p3Ntˆ3Ntq P λ O pN λ ˆ3Ntq fi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi fl so that χ " Px " » - - - - - - - - - - - χ ν " P ν x ax x tr x m x h χ λ " P λ x ax fi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi fl " » - - - - - - - - - - - χ θ χ λ fi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi fl , (3.21)
where P ν (of size p2N t ´Nλ q ˆ2N t ) extracts the unconstrained states, P λ (of size N λ ˆ2N t ) will give the constrained (algebraic) states, x tr , x m and x h are left unchanged, χ ν denotes the new axial inertial components (vector of size N t ´Nλ ) and χ λ is the vector of the constrained states (of size N λ ).

With this new parametrisation and given the form of the coordinate change, the Hamiltonian reads H 1 pχq " H ax pχ ν , χ λ , x m q `Htr px tr , x m q `Hcomp px comp q (3.22)

where H tr and H comp remain unchanged and H ax is H ax pχ ν , χ λ , x m q " H ax px ax "

« P ν P λ ff ´1 » - χ ν χ λ fi fl , x m q " 1 2 » - χ ν χ λ fi fl « P ν P λ ff ´ Q ax px m q « P ν P λ ff ´1 Qχpxmq » - χ ν χ λ fi fl
where we used the matricial notations of Eq. (3.18) and Eq. (3.21).

Q χ px m q (of size 2N t ˆ2N t ) is divided according to the size of χ θ and χ λ such as

Q χ px m q " « Q νν px m q Q νλ px m q Q νλ px m q Q λλ px m q ff
where Q νν (resp. Q λλ ) is a p2N t ´Nλ q (resp. N λ ) square matrix and Q νλ is of size p2N t ´Nλ ˆNλ q.

The constraint equation reads e λ pχ ν , χ λ q " Q νλ px m qχ ν `Qλλ px m qχ λ " 0 ùñ χ λ " ´Q´1 λλ px m qQ νλ px m qχ ν , (3.23) if Q λλ px m q is invertible. For the two-tract example, this simplifies to a single equation which can be solved by hand, thus regenerating the equivalent component method (as in Ref. [START_REF] Wetzel | Power balanced time-varying lumped parameter model of a vocal tract: modelling and simulation[END_REF]).

Using the matricial formulation and taking advantage of the structure of the coordinate change, the implicit set of constraint equations is simplified to a linear set of equations with respect to χ ν . The semi-explicit DA-pH formulation of the full vocal tract is

» - - - - - F VTθ W VT 0 ´y fi ffi ffi ffi ffi fl " » - - - - - - J θ VT G θ w,VT 0 G θ p,VT ´Gθ w,VT 0 0 0 0 0 I pN λ ˆNλ q 0 ´Gθ p,VT 0 0 0 fi ffi ffi ffi ffi ffi fl ¨» - - - - - E VTθ Z VT pW VT q e λ u fi ffi ffi ffi ffi fl (3.24)
where F VTθ " r 9 χ ν , 9 x tr , 9 x m , 9 x h s . The interconnection matrices J θ VT , G θ w,VT and G θ p,VT are defined in the same manner as in the third step of Sec. 3.2.3. The third line is equivalent to Eq. (3.23).

We thus have shown how the use of a permutation matrix simplifies the structure of the coordinate change structure (Eq. (3.21)), and how the Hamiltonian is changed by the coupling procedure (Eq. (3.22) and below). This enabled us to reduce the constraint equation (3.24) to the smallest size possible.

Conclusion

In this section, we assembled a full vocal tract by interconnecting a collection of tracts, the causality of which led to incompatible connections. To address this issue, we proposed a general (automated) assembly method. First, a collection of pHs, the assembly configuration of which is described by a directed graph, has been formulated as a DA-pH. Then, using a method of the litterature, the (implicit) DA-pH has been cast as a semi-explicit dynamical system to allow for simulations.

Then, we applied this procedure to the vocal tract, and took advantage of the particularity of the tract model (namely the canonical form and the matricial expression of the Hamiltonian) to simplify the procedure and reduce the set of implicit constraint equations into a simple matrix inversion problem (linear w.r.t. the constrained DoF).

summary and perspectives

In this chapter, we proposed a new power-balanced FSI vocal tract model, composed of a collection of elementary tracts interconnected by algebraic equations (constraints). Built on the principles of fluid mechanics and thermodynamics, the tract model remains simple and physically interpretable. This was demonstrated using equivalent electrical circuits, which highlighted the existence of a coupling component (gyrator) linking the transverse and axial motion of the fluid. To adopt a nonlinear acoustical description, we proposed an equivalent pHs, the physical variables of which are written as fluctuations around atmospherical conditions (without making any additional assumptions such as a linearising hypothesis). The assembly of the vocal tract model was carried using a new assembly method that relies on directed graphs and DA-pHs.

The model we established allows for a unified representation of the vocal apparatus, including the larynx. Note that applying an incompressible flow assumption (common for the modelling of the glottal flow) regenerates the model of Lopes and Hélie [START_REF] Lopes | Energy Balanced Model of a Jet Interacting With a Brass Player's Lip[END_REF] and Hélie and Silva [START_REF] Hélie | Self-oscillations of a vocal apparatus: A porthamiltonian formulation[END_REF]. Including a passive and simple noise source model and a simple loss model would improve the behaviour of the flow. The considered cartesian geometry is the simplest choice, but it lacks physical realism, whereas building an equivalent tract model in the polar coordinate referential could yield a more realistic description of the flow, especially when a tract narrows. But, as demonstrated in Ref. [START_REF] Arnela | Influence of Vocal Tract Geometry Simplifications on the Numerical Simulation of Vowel Sounds[END_REF], small discrepancies between both models should be observed regarding the acoustical behaviour. An alternative is to consider elliptical shaped tracts. Using a more complex mechanical model of the muscles and tissues surrounding the cavities (see Ref. [START_REF] Alejandro | Port-Hamiltonian Modelling of Fluid-structure Interactions in a Longitudinal Domain[END_REF]) would improve the realism of the vocal tract model, and allow for the use of existing articulatory models (like in Ref. [START_REF] Birkholz | Construction And Control Of A Three-Dimensional Vocal Tract Model[END_REF]).

An interesting approach would be to identify the compatible and incompatible connections within the directed graph, and to automatically derive the constraint matrix C λ (associated with incompatible connections) and interconnection matrix S (associated with compatible connections). The proposed assembly method does not discriminate between the linear and nonlinear constraints. Separating the linear constraints from the nonlinear constraints would allow for the explicit and automatic solving of the former, thus reducing the size of the constraint system of equations. For this last perspective, we point to the work of Beattie et al. [START_REF] Beattie | Port-Hamiltonian descriptor systems[END_REF] and Gernandt, Haller, and Reis [START_REF] Gernandt | A Linear Relation Approach to Port-Hamiltonian Differential-Algebraic Equations[END_REF].

N U M E R I C A L E X P E R I M E N T S

In this chapter, we present numerical experiments to examine the capacities and limitations of the vocal tract model introduced in the previous section. These simulations are carried out using the improved discrete gradient (IDG) method (see chapter 2) and implemented using the Python library PyPHS [Fal].

As PyPHS does not currently support the solving of constraints, we modified the automatically generated C++ simulation code so that, at each time step, the constraint equations are solved prior to the dynamics of the system. The numerical method unfolds as follows:

1. solve Eq. (3.23) for χ λ ; 2. apply the discrete gradient method (see Sec. 2.3.3) on the remaining ODEs.

This simple two-step method enables us to perform simulations for a full vocal tract.

The chapter is organised as follows. In section 4.1, we introduce a few pHs components (structure, fluid losses and radiation condition) to complete the vocal tract, and to enable us to control the geometry. In section 4.2, we examine at a toy model, a simple two-tract system and consider a dynamical setting. In section 4.3, we examine the impact of the number of tracts on the acoustical resonances of a uniform duct. In section 4.4, we consider a simple coarticulation that leads, in section 4.5, to the synthesis of a coarticulation, using a excitation signal.

completing the vocal tract

We introduce a simple mechanical model on the walls, a radiation condition compatible with fluid mechanics, and a fluid loss component. In this chapter, we consider the shifted version of the tract model introduced in Sec. 3.1.5.

Mechanical Model of the Walls

At the wall, the tract model introduced in chapter 3 is driven by an input force. Given this input configuration, we are unable to control the height of each tract. Moreoever, the tract model model does not account for the mechanical behaviour of the structure (i.e. the tissues surronding the vocal tract). Therefore, we connect a mechanical component to the tract wall, which also enables us to drive the walls with an input velocity. We detail this component below.

We consider the damped elastic behaviour of the structure under a low frequency assumption, and neglect its inertia. This leads to the spring-dashpot system depicted on Fig. 23. The state vector x mec is the elongation of the spring ξ around its rest state. The Hamiltonian is H mec pξq " 1 2 kξ 2 where k is the stiffness of the spring and its effort reads e ξ " kξ " F k . The dashpot is characterised by its velocity f d " 9 ξ and the resulting force is af d where a is a dissipation coefficient. Typical values (for e.g. the relaxed cheecks) are k " 845 N ¨m´1 and a " 0.8 kg ¨s´1 [START_REF] Ishizaka | Direct determination of vocal tract wall impedance[END_REF]. For each tract, a mechanical system is connected: to the fluid by the force F ext and velocity v ext ; and to the muscle (the articulator) by the force F m and the velocity v m (see Fig. 23). The velocity of elongation of the spring is 9 ξ and equals the difference of velocity between v ext and v m such as 9

k 1 r 1 v m F m v ext F ext
ξ " ´vext ´vm . The velocity of the dashpot equals the one of the spring. The output force at the interface with the fluid reads F ext " e d `ek . Conversely, for the other interface, we have F m " e d `ek . The full port-Hamiltonian formulation reads

x mec " ξ, H mec pξq " 1 2 kξ 2 , e ξ " kξ and ¨9 ξ w F ext F m ‹ ‹ ‹ ‹ ' " » - - - - - - 0 0 ´1 ´1 0 0 ´1 ´1 `1 `1 0 0 `1 `1 0 0 fi ffi ffi ffi ffi ffi fl ¨eξ " kξ zpwq " aw v ext v m ‹ ‹ ‹ ‹ ' .
The height of a tract can now be controlled through the velocity v m .

Radiation Condition

We now add a radiation condition at the end of the vocal tract. As the tract model accounts for fluid mechanical phenomena, we cannot use standard acoustic impedance types of models. To model the matter exiting the vocal tract, we add a memoryless component, adapted from Ref. [START_REF] Hélie | Self-oscillations of a vocal apparatus: A porthamiltonian formulation[END_REF], that dissipates the kinetic energy of a jet. This component is described by w turb and the non linear law zpw turb q such that w turb " qRk , zpw turb q " 1 2ρ

ˆwturb p hk `h0k qL 0k ˙2 Θpw turb q (4.1)
where Θp¨q is the step function, q Rk the mass flow exiting the vocal tract, and k the index of the tract, to the right interface of which we connect this lossy component (lips, nose). Note that the power dissipated by the component is positive or zero as P turb " w turb zpw turb q " w turb

1 2ρ ´wturb h k L 0k ¯2 Θpwq ě 0.

Losses of a Laminar Flow in a Rectangular Duct

Following the approach of Maeda [START_REF] Maeda | A digital simulation method of the vocal-tract system[END_REF], we add a simple fluid loss component. Without this component, and under a glottal excitation, the axial velocity in the vocal tract would eventually reach infinity. As in Ref. [START_REF] Maeda | A digital simulation method of the vocal-tract system[END_REF], we consider the formula of the laminar resistance in a rectangular duct given in Ref. [START_REF] Kenneth | Airflow and Turbulence Noise for Fricative and Stop Consonants: Static Considerations[END_REF] R µ p hk q " 8µ 0 S k p hk q (4.2)

where R µ p hk q is the resistance per unit length, µ 0 is the dynamic viscosity of the air, and Sp hk q " p hk `h0k qL 0k is the cross section of the k-th tract.

We describe this fluid loss formula as a memoryless component of variable w pkq µ " qk and law z µ pw pkq µ q " R µ p hk q 0k w pkq µ , (4.3) the associated power P µ " z µ pw pkq µ qw pkq µ " R µ p hk qpw pkq µ q 2 of which is always positive or null.

This component is only used in section 4.3, section 4.4 and section 4.5. The equivalent electrical circuit of the shifted tract model, augmented with the resistors corresponding to the dissipative component introduced above, is shown on Fig. 24 m + m Πy Connecting a spring-damper mechanical system the external port of the wall pF W , v W ) is equivalent to connecting a capacitor and a resistor in series with the fluid mass Πy (at the top of Fig. 24).

v W Πy m + m - ṁ Πy m + m ∂ m H = -ψtot ṁ R µ ( h) z µ (w µ ) w µ νL q L νR q R R µ ( h) z µ (w µ ) w µ ∂ h H = S
We may now proceed with performing numerical experiments.

the two-tract system: a toy model

To examine the model and its signal in a somewhat exhaustive way, we start by the most basic configuration: a simple two-tract system of total length 17 cm, as the one depicted on Fig. 25 (where we omitted the " notation for legibility). We connect two spring-dashpot mechanical systems to the walls and a radiation component to the rightmost surface of T 2 (S L2 ). We examine the behaviour of the model for a simple dynamical configuration, where the wall of T 2 is driven by an input velocity, and focus on: the accumulation of mass due to the change of geometry; the power-balance; and the existence of acoustical resonances.

q R 1 = q L 2 T 1 T 2 {h 1 , m 1 , Π y1 } {h 2 , m 2 , Π y2 } v W 1 e y F W 1 e y v W 2 e y F W 2 e y ψ L1 q L1 e x ψ R2 q R2 e x
The simulation is split into three temporal phases:

(1) a static phase t P rt 0 , t 1 " 1 ss where both tracts are identical in geometry (h 1 " h 2 " 1cm);

(2) a dynamical phase t P rt 1 , t 2 " 1.05 ss where the wall of tract T 2 is driven to reach a target height h 2 " 2 cm;

(3) a final static phase t P rt 2 , t max " 3 ss similar to phase (1).

To produce frequency representations, an enthaply impulse is fed to the vocal tract through the leftmost surface of T 1 (S L1 ) at t imp 1 " 0.3 s for phase (1), and at t imp 2 " 2 s for phase (3). In App. C (see page 131), we plot the velocity input signal used to drive the wall of T 2 (see Fig. 45) and the enthalpy impulse (see Fig. 46). The physical and numerical parameters are detailed in Tab. 3 and Tab. 4. The physical parameters of the mechanical systems connected to the walls will remain identical for every simulation in this chapter. The value of the stiffness of the spring was taken from Ref. Ishizaka, French, and Flanagan [START_REF] Ishizaka | Direct determination of vocal tract wall impedance[END_REF] (relaxed cheeks setting). On all plots, t 1 and t 2 are indicated with vertical red dashed lines. On the relevant figures, the static phases (1) and (3) are denoted with a blue background, and the dynamical phase (2) is denoted with an orange background.

We first focus here on the dynamical phase (2) of the simulation. On Fig. 26, we plot the heights h 1 (plain green curve) and h 2 (dashed black curve) of both tracts, and the total mass m 1 `m2 versus time during the movement. As tract T 2 becomes larger, mass is being drawn from the exterior into the system, thus increasing the total amount of mass in the vocal tract. Now, we assert the passivity of the simulation. On Fig. 27, we plot the local error on the power-balance for the full simulation (at the top) and only for the dynamical phase (2) (at the bottom). The local error on the power balance remains within the precision of the chosen float-point representation and the numerical conditioning of computations. The error is maximal during the dynamical phase, where the system exhibits the most nonlinear behaviour (mass convection).

t

Local error (W)

Figure 27: Local error on the powerbalance for the two-tract model. At the top, t P r0, t max s.

At the bottom, we only show the dynamical phase (2) with t P rt 1 , t 2 s.

We now consider the frequency responses of the static configurations associated with phases (1) and (3). Here, the frequency response (see remark 8) is the ratio of the Fourier transforms of the flowrate u lips pf q " q L1{ρ 1 and u g pf q " q R2{ρ 2 respectively at the leftmost surface S L1 (the glottis) and rightmost surface S R2 (the lips). The peaks of the frequency response correspond to acoustical resonances, called formants for time-invariant assumptions. We plot the frequency responses for the static phases (1) and (3) on Fig. 28. Frequency response H(f ) = u lips /u g

Step (1)

Step (3) At phase (1), the system shows a main resonance at f 1 " 130Hz. In phase (3), this formant shifts to f 1 " 168 Hz as a consequence of the wall movement in phase (2). The amplitude of the resonance and the (pseudo) quality factor diminish because of the height expansion of T 2 , which increases the magnitude of the radiation effect.

For both static configurations, the frequency respresentations show a secondary resonance, the magnitude of which is small compared to the one of f 1 ( f 2 « 400 Hz for phase (1) and f 2 « 500 Hz for phase (2)).

The presence of two resonances indicate that a given vocal tract, made of N t , will exhibit N t formants. This behaviour is investigated in the next section.

Remark 8 (frequency representation): Both plots on Fig. 28 are not transfer functions, because the vocal tract model is nonlinear. However, as we are looking at small fluctuations around static configurations (phases (1) and (3)), the computed frequency responses are close to the transfer functions of the linearised model. As expected, the numerical vocal tract model accounts for variations of geometry, and mass, and exhibits acoustical resonances while preserving the power balance. However, a straight, uniform, 17-centimeter-long acoustic tube should exhibit a first resonance at around f " 500 Hz. In our case, in phase (1), where the geometry of the vocal tract is the one of a uniform tube, the model exhibits a resonance at f " 130 Hz. This discrepancy is due to the lumped-parameter modelling combined with the low number of tracts that we used to discretise the whole vocal tract. In the next section, we examine the impact of the number of tracts on the acoustical resonances of the model.

impact of the discretisation on the acoustical behavior

We now consider the case of a uniform duct of total length L tot " 17 cm, the acoustical resonances of which are odd multiples of the fundamental frequency f 1 " 500 Hz (for an ideal closed-open resonator with an ideal radiation condition, see Ref. [START_REF] Chaigne | Acoustique des instruments de musique (2e {é}dition revue et augment{é}e)[END_REF]3.5.2,p. 273]). To study the impact of the number of tracts on the frequency location of the formants, we discretise the vocal tract with a range of N t P 2, 22 identical tracts of equal length, height and width, such that 0 " Ltot {2Nt and h 0 " L 0 " 1 cm.

To estimate the frequency response for each discretisation setup, we feed an enthalpy impulse to the leftmost surface S L1 at t " t max {3 with t max " 0.1 s. The sampling frequency is identical for each simulation (5 MHz). The remaining physical and numerical parameters are identical to those of section 4.2.

On Fig. 29, we plot the frequency response for each discretisation setting, from N t " 2 all the way to N t " 22. The number of formants that the vocal tract exhibits increases with the number of tracts, so that each new tract adds a new resonance. As the total number of tracts increases, the frequency location of each formants converges to their corresponding theoretical value (vertical orange dashed lines). The discrepancy between the theoretical values and the numerical results increases with the frequency. Note that the first formant is always underestimated, whereas the higher ones are always overestimated. This is most likely caused by the (nonlinear) lumped radiation condition that we added at the open end of the resonator, acting as a correction length.

The plot also shows that considering 15 tracts to discretise the vocal tract is enough to reconstruct the three first resonances of the acoustical system.

a simple coarticulation

We now consider the case of an articulating vocal tract, the shape of which varies with time. Reflecting on the previous numerical experiment, we use a total of N t " 15 tracts to discretise the resonator. This setting allows us to examine some acoustical and energetical phenomena in a single simulation.

To each tract, we add the simple model of laminar resistance [START_REF] Kenneth | Airflow and Turbulence Noise for Fricative and Stop Consonants: Static Considerations[END_REF] (detailed in Sec. 4.1 and used in Ref. [START_REF] Maeda | A digital simulation method of the vocal-tract system[END_REF]) to yield a more realistic behaviour and more stable simulations. Indeed, the resulting 15-tract systems involves 14 constraints (one per connection), thus increasing the stiffness of the dynamical system. Adding supplementary losses helps the Newton-Raphson solver to converge. The stiffness of the numerical problem requires a sampling frequency of 1 MHz. Moreover, solving the constraints prior to the dynamics, instead of solving the implicit system altogether, is not ideal. This numerical aspect was not investigated further and will be discussed in the conclusion as a perspective.

f 1 = 656Hz f 2 = 1580Hz f 1 = 155Hz f 2 = 2000Hz
Frequency response H(f ) = u lips (f )/u g (f )

[a]

[i] In a real life setting, the expected formant locations for a vowel [a] are around f 1 « 850 Hz and f 2 « 1610 Hz and f 1 « 240 Hz and f 2 « 2400 Hz for vowel [i], according to Ref. [START_REF] Catford | A Practical Introduction to Phonetics[END_REF]. In our case: for vowel [a], f 1 « 656 Hz and f 2 « 1580 Hz; for vowel [i], f 1 « 155 Hz and f 2 « 2000 Hz. The discrepancies are most likely due to the simple geometrical hypotheses that come with a lumped-parameter approach.

The magnitude of the resonances are higher in vowel [a] than for vowel [i]. This is a consequence of the fluid dissipation model we introduced, the dissipation rate of which increases as the cross section of a tract diminishes. In the vowel [i], toward the mouth, the cross sections of the tracts are smaller (see Fig. 30). Therefore, the dissipation rate is higher, leading an output signal with less magnitude than for vowel [a], where the mouth is wide open. This difference in magnitude has an impact on voiced signals, that we examine in the next section.

We now compare the magnitude of the power flows induced by the acoustical impulses, fed to the vocal tract at the static phases (1) and (3), to the power flows induced by the movement of the wall at phase (2). The former are related to acoustic effects, whereas the latter are related to fluid mechanical effects.

We thus look at the power flows (time variation of stored energy) of each tract, the product of the effort and the flux variables. At the top of Fig. 32, we plot the power flows for phases (1) and (2), omitting phase (3) as it is redundant with phase (1). Each colored curve denotes a tract, the dark blue curve being T 1 , the dark red curve being T 15 . During phase (1), the peak value of the power flows is approximately 10 ´8 W. During the movement, in phase (2), the power flow is of unitary value, several orders of magnitude above the the acoustical fluctations of step (1).

At the bottom left of the figure, we zoomed in on the instant where the impulse, of unitary magnitude, is fed to T 1 . The power flow, induced by this impulse, propagates from T 1 (dark blue) to T 15 (dark red). Due to the radiation condition, losses in the mechanical systems and in the fluid, the variation of stored energy decreases. Then, the tract exhibits an oscillatory behaviour at its highest resonant frequency (« 15.5 kHz, see bottom right plot). This artefact is most likely due to the use of an inadequate numerical method for the solving of the constraints of a stiff problem. Again, solving the constraints prior to the dynamics of the system is not ideal and probably causes this high frequency oscillation.

With this numerical experiment, we have shown that the articulated vocal tract exhibits acoustical resonances that are realistic. We examined the power flows related to the acoustical and fluid mechanical phenomena, and found that the latter is a few orders of magnitude above the former. However, the results have to be confirmed by new numerical experiments that use an adequate numerical method to properly handle the constraints.

vowel synthesis

In this section, we excite the vocal tract with a simple source and consider the coarticulation between vowel [a] and [i]. As in the previous section, we consider three phases:

(1) a static phase, t P r0, t 1 " 0.75 ss, where the vocal tract articulates the vowel [a];

(2) a dynamical phase, t P rt 1 , t 2 " 1.85 ss, where the walls are driven to reach the geometry of vowel [i];

(3) a static phase, t P rt 2 , t max " 3 ss, where the vocal tract articulates the vowel [i].

The geometric data for the geometry are extracted from Ref [START_REF] Story | Vocal Tract Area Functions from Magnetic Resonance Imaging[END_REF]. The velocity signal used to drive the walls is shown in App. C (see Fig 48).

We consider a simple excitation signal constituted of replicated elementary impulses as the one used in the previous section (see Fig. 46 in App. C), so as to obtain a train of impulsion at base frequency f 0 " 200 Hz.

We plot the simulated flowrate at the lips and its time-derivative on Fig. 33. Flow rate at the lips and its time derivative The original flowrate exhibits a prominent low-frequency component (top of Fig. 33), due to the phenomena of mass convection. To isolate the acoustic components, we plot the time derivative of the flow rate (see bottom plot). The amplitude of the signal is greater for the [a] consonant, as the mouth is open, and smaller for the [i] consonant. This discrepancy is mainly due to the loss model and to the difference of surface area of the open mouth.

On the time derivative of the flowrate, the low frequency transients are still present and can be seen on the spectrogram of the signal, shown on Fig. 34 (see the low frequency range juste after t 1 and jud after t 2 ). For this figure, we describe each plot from top to bottom. At the top, we plot the geometry of the vocal tract at the start and at the end of the simulation. Below, we plot an articulation coefficient, the average normalised distance 1 between the starting and ending geometric configurations: this coefficient is null in phase (1), is equal to one when in phase (3), and the articulation coefficient varies linearily with time in phase (2). Note that it could be possible to use a more realistic signal, extracted from resonant magnetic imaging for instance. Below, we plot the total transverse flow energy, the energy-normalised spectrogram of the output signal (time derivative of the flowrate at the lips), and the energy per time frame of the original spectrogram (used to normalise the spectrogram).

On the spectrogram, the formants are clearily visible and move with the same timing as for the articulation coefficient, indicating a correlation between the articulation and the frequency shift of the acoustical resonances. The low frequency transients can be seen at the bottom of the spectrogram, at t " t 1 and t " t 2 . As mentionned earlier, the magnitude of the signal is smaller for the vowel [i], as the fluid losses in the mouth are more prominent (see bottom plot of Fig. 34). As expected, the energy of the transverse flow is maximal during the movement. To show the acoustical propagation phenomena, we plot on Fig. 35 the input enthalpy ψL1 at the glottis and the mass rates q1 and q15 at both extremities of the vocal tract. We zoomed in on two impulses. At each impulse, there is a delay between the output mass rate at the glottis and at the lips. The value of this time delay correspond approximately to the time it takes for the enthalpy impulse to travel along the vocal tract (« 5 ms).

1 1 {Nt ř Nt i"1 |hiptq´hipt0qq {|hpt0q´hiptmaxq|

summary and perspectives

In this section, we completed the vocal tract model with additional components and performed numerical experiments of increasing complexity. Using a simple twotract model, we have shown that the numerical vocal tract model handles mass accumulation and geometry variations, and exhibits acoustical resonances. Then, we investigated the impact of the number of tracts on the formants by discretising a uniform duct with an increasing amount of tracts. The frequency location of each formant converges to the expected value as the number of tracts increases. We then studied a practical coarticulation example, and proposed an example of synthesis.

The above numerical experiments are not easy to perform, as the numerical method is not well suited to handle constraints. To address this issue, the constraints and dynamics have to be simultaneously solved by the Newton-Raphson solver. This would allow for the use of a smaller sampling frequency. We refer to the work of Müller [START_REF] Müller | Time-continuous Power-balanced Simulation of Nonlinear Audio Circuits : Realtime Processing Framework and Aliasing Rejection[END_REF] for a class of power-balanced numerical methods that would enable us to use a lower-sampling frequency.

Moreover, the time needed for the computations is too large to consider a real-time implementation. For a system composed of 15 tracts, a simulation of 3 seconds has to run for an hour. This is most likely due to: the C++ code generated by PyPHS that is not optimised; and the very high sampling frequency (1 MHz) needed to solve the stiff problem. However, the tract model may be compatible with a quadratisation method [START_REF] Lopes | Approche passive pour la modélisation, la simulation et l'étude d'un banc de test robotisé pour les instruments de type cuivre[END_REF], which would allow for an explicit formulation of the numerical problem. The presence of additional constraints has to be adressed by extending the quadratisation method to the case where the Hamiltonian is parametrised by instantaneous constrained variables. Finally, the above model can be used with an articulatory model to perform more complex coarticulation tasks. A more thorough work has to be conducted to fully understand the capacities of our vocal tract model, and identify the possible optimisation that can be applied to it. In this chapter, we are interested in the interconnection of three resonators with a small common volume called a three port junction, shown on Fig. 36. At this intersection, the cavity is submitted to the pressure and flow rates of each resonator. This is the case in the vocal tract, where the lower airways split at the nasal branching into the mouth and nasal cavities, and in woodwinds instruments, where lateral holes can be considered to be short open resonators connected all along the main bore. Several models exist for such a junction, some of them are mentioned below.

A PA S S I V E T H R E E -P O RT J U N C T I O N F O R F L U I D M E C H A N I C S A N D A C O U S T I C S

The most basic junction model (used in chapter 3) assumes the continuity of pressure forces and that the sum of incoming mass fluxes at the junction vanishes. However, this basic model does not account for any geometrical parameters. This issue has been adressed in the linear acoustical domain in Refs. [START_REF] Vincent Dubos | Theory of Sound Propagation in a Duct with a Branched Tube Using Modal Decomposition[END_REF] and [START_REF] Chaigne | Acoustique des instruments de musique (2e {é}dition revue et augment{é}e)[END_REF]Sec. 7.7]: a model of a branched resonator is derived using Green functions, a variational formulation and a modal decomposition for the 2D rectangular geometry1 , resulting in a purely inertial model. The parameters of the model are function of the geometry of the junction and the mass of fluid connecting volume. This model has no guarantee w.r.t. the passivity, due to the presence of negative masses, the value of each is derived numerically. In the fluid mechanical domain, Mora Araque [START_REF] Alejandro | Port-Hamiltonian Modelling of Fluid-structure Interactions in a Longitudinal Domain[END_REF] proposed a purely compressible node model that connects two incompressible pHs models. This is, to our knowledge, the only model of junction within the pH framework. The geometrical parameters of the model are small and delicate to estimate, while their importance in the numerical simulations is significant.

In summary, the first model does not guarantee the passivity, and the second is not suited for the interconnection of acoustical components and accounts for only two connections. As there is a growing interest for physical guaranteed-passive modelling of linear and nonlinear acoustical systems [BH16; TC20; HR16] and fluid mechanical systems [CML19; HS17; Hau+19; Mor20], the need for a passive and physically interpretable junction model that can handle all these different phenomena increases.

To this end, we propose a simple and passive three-port junction model composed of a 2D parallelepiped volume of homogeneous mass density. It can be connected to other dynamical systems through three open surfaces. We start by establishing a macroscopic base model (M1) using a set of geometrical and physical hypotheses. From this fluid mechanical model, we derive five different pHs that account for (M1ab) compressible fluid mechanics, (M1c) fully and (M2) weakly nonlinear acoustics, and, at last, (M3) linear acoustics (see Fig. 37). (M2) and (M3) incorporate linearising approximations that neglect, respectively, the nonlinearities of the state law and of the mass convection.

M1

M1a

M1b M1c

M2 M3

Fluid mechanics Nonlinear acoustics Linear acoustics Throughout this iterative process, we give some physical interpretations and analogies with classical approaches through the use of equivalent electrical circuits and state-space representation. The passivity of every model is ensured by the use of pHs.

We apply the same modelling methodology as in chapter 3, to which we add elements regarding the generalisation of our approach. This would allow for the use of more complex approximating basis functions. The chapter is organised as follows. In section 5.1, we describe the considered geometry and remind the reader of the physical hypotheses that we consider. In section 5.2, we discretise the model in order to obtain a finite dimensional model. Then, in section 5.3, we cast the model into the pHs framework and derive 5 different models ranging from fluid mechanics to linear acoustics. Finally, in section 5.4, we compare models together and give some physical interpretations of the linear model by examining its state-space representation.

hypotheses and approximations

The three-port junction of Fig. 36 is modeled as a (H0) 2D rectangular time-invariant volume, depicted on Fig. 38, delimited by the spatial domain Ω " px, y, zq P r´ 0 , 0 s ˆr0, h 0 s ˆr0, L 0 s

( of boundary BΩ " S L Y S b Y S R Y S T .
The surface S b (located at y " 0) is an impervious boundary while S L , S R and S T (respectively located at x " ´ 0 , x " 0 and y " h 0 ) are open boundaries. The volume V 0 of the junction is given by V 0 " 2 0 L 0 h 0 .

q T ψ T ψ L q L ψ R q R -0 0 x = 0 0 h 0 S R S L S T S b e y e x

Ω

Figure 38: The considered 2D cartesian geometry. To be compared with general case of Fig. 36 where: system 1 is connected to S L , system 2 to S R and system 3 to S T .

We remind the reader of the physical hypotheses we consider (see chapter 1):

(H1) the fluid is assumed to be perfect (no viscosity); (H2) the flow is irrotational; (H3) the process is isentropic and adiabatic;

(H4) the heat capacity at constant volume C V is constant; (H5) the gas is assumed to satisfy the perfect gas assumption;

to which we add the (H6) lumped-parameter hypothesis:

(H6) the volumetric mass density ρ is homogeneous within Ω.

Hypotheses (H1-5) describe a lossless fluid mechanical system, the dynamics of which is governed by the conservation of mass and momentum Bρ Bt " ´divpρvq and Bv Bt " ´gradpψq the latter being shown in its specific form2 and where ψ " 1 2 v ¨v ` pρq `ρ B pρq Bρ is the specific enthalpy and pρq is the specific internal energy of the fluid.

We now detail the boundary conditions. For each S i , i P tL, R, T u, we define the surface averaged enthalpy ă ψ ą S i , and its power-conjugated quantity, the mass flow

q S i ă ψ ą S i " 1 {S i ij S i ˆ1 2 v ¨v ` pρq `ρ B pρq Bρ ˙dS and q S i " ij S i pρvq ¨n dS, (5.1)
where n is the outward unit vector, orthogonal to the associated surface. The total energy E stored by the fluid in Ω is the sum of a kinetic and thermodynamic (or compression) term such as

Epv, ρq " ¡ Ω ˆ1 2 ρv ¨v `ρ pρq ˙dV , (5.2)
where, under (H3-5), the internal specific energy pρq reads (see see Eq. (B.1.2)3 in App. B.1)

pρq " P 0 pγ ´1qρ 0 « ˆρ ρ 0 ˙γ´1 ´γff .
(5.3)

All the above equations and hypotheses constitutes the physical description of the junction. We now establish the first macroscopic finite dimensional model.

spatial discretisation: deriving the macroscopic model (m1)

The macroscopic model (M1) is obtained by discretising the system similarily as in the FEM method [START_REF] Ern | Theory and Practice of Finite Elements[END_REF]. We project the velocity v and volumetric mass density ρ on a basis function matrix P of size p3 ˆnq which performs the following approximations

« v ρ ff « P « v ρff " « P v O p2ˆnρq O p1ˆnvq P ρ ff « v ρff
where v " vptq and ρ " ρptq are finite dimensional vectors respectively of size n v and n ρ with n " n v `nρ . They contain the degrees of freedom used to approximate each infinite dimensional physical variable. P is divided accordingly to the dimension of v, ρ, v and ρ so that P v is of size p2 ˆnv q and P ρ is of size p1 ˆnρ q. Each component of P is a shape function.

Substituing v and ρ by their discretised counterparts in the energy function of Eq. (5.2) yields Epv, ρq « Êp v, ρq "

¡ Ω "
1 2 pP ρ ρq p v P v q pP v vq `pP ρ ρq pP ρ ρq  dV .

(5.4)

The weak form of the conservation of mass and momentum equations is (see App. B.2 for details)

d dt ¡ Ω ϕρ dV " ´ij BΩ pϕρvq ¨n dS `¡ Ω gradpϕq ρv dV , (5.5a) d dt ¡
Ω ϕv dV " ´¡ Ω divpϕψq dV `¡ Ω divpϕqψ dV , (5.5b)
where ϕ : R 2 Ñ R, with ϕpx, yq P C 8 , is a smooth function w.r.t. x and y.

As in chapter 3, we make the simplest choice of basis that is compatible with our hypotheses and boundary conditions.

Definition 1 (Space of polynomials): We denote by P k with k P N the space of polynomials of degree less or equal to k.

We choose the space of affine functions P 1 for the velocity and (H6) piecewise constant P 0 for the mass density, and build P as where, as in chapter 3, v mx ¨ex accounts for a uniform axial motion, vcx { 0 ¨ex for a velocity gradient and vyy {h 0 ¨ey for the transverse flow. In the following, we omit the time dependency of these quantities.

Remark 9 (Transverse velocity): The transverse velocity is only discretised with the single shape function y {h 0 as the velocity of the fluid must be zero at the bottom wall (no uniform transverse motion of the fluid).

Then, the energy is given by

Êp v, ρq " 1 2 v Qp ρq v `ρ p ρqV 0 where Qp ρq " ρ ¡ Ω pP v P v q dV " ρ ¡ Ω ¨1 ´x 0 0 ´x 0 x 2 2 0 0 0 0 y 2 h 2 0 ‹ ‹ ‹ ' dV " ¨V0 ρ 0 0 0 V 0 ρ 3 0 0 0 V 0 ρ 3 ‹ ‹ ' ,
and p ρq is given by Eq. (5.3).

Remark 10 (Higher order of discretisation): The above computations can be performed with higher order of basis function providing a suitable form function matrix.

We now derive the dynamics equations that governs each component of v and ρ. In Eq. (5.5), we replace the variables v and ρ by their approximation P v v and P ρ ρ " ρ and replace ϕ by one of the shape function of P in Eq. (5.6). For instance, choosing ϕ " 1, the shape function associated with ρ, in the conservation of mass yields d dt where we see that only the affine component v c and v y account for the compressibility. As v mx accounts for an incompressible axial flow, it does not appear in this equation, as it is related to the compressible behaviour of the fluid.

¡ Ω ρ dV " ´ij BΩ p ρP v vq
For the components v, we now detail the computations for the first component v mx , associated with the shape function 1 e x . Setting ϕ " 1 e x in Eq. (5.5b) gives d dt

¡ Ω P v v e x dV " ´ij BΩ ψ n dS
where ψ is the discretised enthalpy obtained by replacing v and ρ by their approximation in Eq. (5.1). Computing the volume integral on the left hand-side and expanding the surface integral on the right hand-side yields

d dt pV 0 v mx q " ´ij S L ψ dS `ij S R ψ dS " L 0 h 0 `ă ψ ą S R ´ă ψ ą S L
where the surface-averaged enthalpies are unknown inputs and S " L 0 h 0 denotes the cross section of the junction.

For v c and v y , we respectively choose the shape functions ϕ " t ´x 0 e x , y h 0 e y u which gives

d dt p V 0 3 v c q " S ă ψ ą S L `S ă ψ ą S R ´V0 2 0 ˜2p ρ B p ρq B ρ ` p ρqq `v2 c 3 `v2 mx `v2 y 3 ḑ dt p V 0 3 v y q " S ă ψ ą S T `V0 2 0 ˜2p ρ B p ρq B ρ ` p ρqq `v2 c 3 `v2 mx `v2 y 3 ¸.
At last, we approximate each output mass flow of Eq. (5.1) according to Eq. (5.6) such as q S L « qL " ´L0 h 0 ρ pv c `vmx q , q S R « qR " L 0 h 0 ρ p´v c `vmx q and q S T « qT " 2L 0 0 ρv y .

(5.9)

We can now cast this macroscopic model as a pHs.

ph formulations: from (m1a) compressible fluid mechanics to (m3) linear acoustics

In this section, we provide five different pH models of the three-port junction model. (M1a-c) are equivalent: only their formulation differs. (M2) is equipped with a linear state law, and (M3) is fully linear.

Model (M1a): Natural Parametrisation

We first choose v and ρ as energy variables, as they parametrise the chosen kinematic and thermodynamical spatial approximations (see Eq. (5.6)). Then, the state vector is x 1 " rv mx , v c , v y , ρs and the associated Hamiltonian is

H 1 px 1 " " v ρı q " 1 2 v ¨V0 ρ 0 0 0 V 0 ρ 3 0 0 0 V 0 ρ 3 ‹ ‹ ' v `V0 P 0 γ ´1 "ˆρ ρ 0 ˙γ ´γ ρ ρ 0  (5.10) leading to the effort vector e 1 , e 1 " ∇ x 1 H 1 px 1 q " » - - - - - e vmx e vc
e vy e ρ fi ffi ffi ffi ffi fl

" » - - - - - - - ρV 0 v mx ρV 0 vc 3 ρV 0 vy 3 1 2 V 0 ´v2 mx `v2 c 3 `v2 y 3 ¯`V 0 P 0 γ pγ´1qρ 0 " ´ρ ρ 0 ¯γ´1 ´1 fi ffi ffi ffi ffi ffi ffi fl
.

Using the equality e ρ " ´ţΩ ψ dV " ´V0 ă ψ ą Ω , we express the outputs and equations of dynamics in terms of the components of e and the inputs ă ψ ą S L , ă ψ ą S R and ă ψ ą S T (App. D). Then, we identify the following algebraic-differential formulation by casting the equation in a matrical form such as

dx 1 dt y 1 ¨dvmx dt dvc dt dvy dt d ρ dt qL qR qT ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' " » - - - - - - - - - - - - - - 1 {2 0 ´1{2 0 ´3{ 0 V 0 3 {2 0 3 {2 0 3 {h 0 V 0 ´3{h 0 3 { 0 V 0 ´3{h 0 V 0 ´1{2 0 ´3{2 0 1 {2 0 ´3{2 0 3 {h 0 fi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi fl S 1 » - - - - - - - - - - - - - - e vmx e vc e vy e ρ ă ψ ą S L ă ψ ą S R ă ψ ą S T fi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi fl e 1 u 1 .
This choice of energy variables renders the interconnection matrix not canonical as S 1 depends on the physical parameters 0 , L 0 and V 0 . In the case of a canonical matrix, the power flows are encoded only with ˘1, which simplifies the interpretation of the interconnection of efforts and flows.

Model (M1b): Canonical Formulation

A canonical formulation of a pHs enables one to write the power exchanges between the components with ˘1 in the interconnection matrix S. It encodes an elementary connection, like the rigid connection between the mass and the spring in the mass-spring-damper example of chapter 2. This property of canonicity relies mainly on the choice of state vector and the presence of additional algebraic components (transformers and gyrators) possibly encoding physical coupling (like in chapter 3).

Here, the canonical formulation is obtained by applying a simple change of coordinates. We define a new state vector x 2 such as

x 2 " Mx 1 " » - - - - - ´ 0 ´ 0{3 0 0 0 ´ 0{3 0 0 0 0 h 0{3 0 0 0 0 V 0 fi ffi ffi ffi ffi fl ¨vmx v c v y ρ ‹ ‹ ‹ ‹ ' ": ¨νL ν R ν T m ‹ ‹ ‹ ‹ ' " ¨´ 0 pv mx `vc {3q 0 pv mx ´vc {3q vyh 0{3 ρV 0 ‹ ‹ ‹ ‹ '
(5.11) so the Hamiltonian reads

H 2 px 2 q " m 2 2 0 pν 2 L `ν2 R `νL ν R q `3m 2h 2 0 ν 2 T `V0 P 0 γ ´1 "ˆm m 0 ˙γ ´γ m m 0  where m 0 " ρ 0 V 0 . The associated effort vector e 2 " ∇ x 2 H 2 px 2 q is e 2 " ¨m 2 0 pν L `νR{2q m 2 0 p ν L{2 `νR q 3mν T{h 2 0 1 2 2 0 pν 2 L `ν2 R `νL ν R q `3 2h 2 0 ν 2 T `P0 γ ρ 0 pγ´1q ´pm {m 0 q γ´1 ´1¯‹ ‹ ‹ ‹ ‹ ‹ ' " ¨q L qR qT ă ψ ą Ω ‹ ‹ ‹ ‹ '
where ă ψ ą Ω is the volume averaged enthalpy.

Substituing x 1 by M ´1x 2 in the algebraic formulation of (M1a) yields the following pH canonical formulation

¨dν L dt dν R dt dν T dt dm dt qL qR qT ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' " » - - - - - - - - - - - - - - - 1 ´1 1 ´1 1 ´1 ´1 ´1 ´1 1 1 1 fi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi fl S 2 » - - - - - - - - - - - - - - - e ν L e ν R e ν T e m ă ψ ą S L ă ψ ą S R ă ψ ą S T fi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi fl
(5.12)

where the exchange of power are indicated by ˘1. This simple structure enables us to build an equivalent electrical circuit, using the fluid mechanical electrical analogy presented in chapter 3. We show the circuit on Fig. 39, where the topology of Fig. 38 has been kept. e m ṁ νL Inductors represent fluid masses and the capacitor represents the compressibility of the fluid. On the differential algebraic formulation, the dynamics of the kinetic components (ν L , ν R and ν T ) are analogous to the Kirchhoff loop rule and the remaining equations are analogous to the Kirchhoff current law [Md 18].

q L νR q R Πy q T ψ L ψ R ψ T
Both left and right branches represent the axial flow (components ν L and ν R ) and the top branch (component ν T ) represents the transverse flow. In the solid mechanical analogy, the mass component m interprets as a "two-dimensional spring" that is submitted to the velocities of the two axial masses and transverse masses.

This model is conservative (it does not exhibits any resistors). The simplicity of the structure is a direct consequence of the canonical interconnection matrix. The dynamics of the system is examined in the last section of the chapter.

For later developpements, we divide the interconnection matrix as

S 2 " « J 2 G p2 ´G p2 O p3ˆ3q ff with J 2 " » - - - - - 0 0 0 1 0 0 0 1 0 0 0 1 ´1 ´1 ´1 0 fi ffi ffi ffi ffi fl and G p2 " » - - - - - ´1 0 0 0 ´1 0 0 0 ´1 0 0 0 fi ffi ffi ffi ffi fl .
(5.13) Both (M1a) and (M1b) are equivalent fluid mechanical models. We now introduce the equivalent model (M1c), the formulation of which is closer to acoustics.

Model (M1c): Nonlinear Acoustics

In this section, we assume that the system is submitted to the atmospheric conditions, characterized by a fluid at rest (no velocity) and at atmospheric pressure p " P 0 . From this assumption, we derive a new model (M1c), the quantities of which (pressure, velocities, volumetric mass density) are expressed as fluctutations around their atmospherical value.

To this end, we apply the shift procedure introduced in chapter 2. Like models (M1a-b), (M1c) is exact, as we do not introduce any linearising approximations.

The set of static inputs that corresponds to the atmospherical conditions is u ‹ " ră ψ‹ ą S L , ă ψ‹ ą S R , ă ψ‹ ą S T s " r0 , 0 , 0s under which we assume that (M1b) is stationary. The associated stationary state vector is x ‹ 2 and verifies 9 x ‹ 2 " 0. Giving u ‹ and x ‹ 2 , the dynamics of ν L , ν R and ν T , on Eq. (5.12)), read e m px ‹ 2 q " 0, and the one of m reads e ν L px ‹ 2 q `eν R px ‹ 2 q `em px ‹ 2 q " 0. For this system of equations, a trivial solution is x ‹ 2 " r0 , 0 , 0 , m 0 s .

From there, we define the shifted energy variable vector is x 3 " x 2 ´x‹ 2 " r νL , νR , νT , ms and the shifted Hamiltonian

H 3 px 3 " r νL , νR , νT , ms q " H 2 px 3 `x‹ 2 q ´e2 px ‹ 2 q ¨x3 ´H2 px ‹ 2 q " m `m0 2 2 0 p ν2 L `ν 2 R `ν L νR q `3p m `m0 q 2h 2 0 ν2 T `P0 V 0 γ ´1 "ˆ1 `m m 0 ˙γ ´γ m m 0 ´1 , Ecp mq
(5.14) leading to the (shifted) effort vector e 3 " ∇ x 3 H 3 px 3 q

e 3 " ¨p m`m 0 q 2 0 p νL `ν R{2q p m`m 0 q 2 0 p νL{2 `ν R q 3p m`m 0 q νT{h 2 0 1 2 2 0 p ν2 L `ν 2 R `ν L νR q `3 2h 2 0 ν2 T `P0 γ ρ 0 pγ´1q ´pp m`m 0 q {p m`m 0 q 0 q γ´1 ´1¯‹ ‹ ‹ ‹ ‹ ‹ ' " ¨q L qR qT ă ψ ą Ω ‹ ‹ ‹ ‹ '
, the three first components of which are identical to the one of (M1b) and ă ψ ą Ω is the fluctuation of enthalpy around its state at rest value.

The algebraic differential formulation is

¨dx 3 dt y 3 '" S 3 ¨e3 u 3 ' (5.15)
with S 3 " S 2 and where u 3 and y 3 are the shifted inputs and ouputs with u 3 " ră ψ ą S L , ă ψ ą S R , ă ψ ą S T s and y 3 " y 2 " r qL , qR , qT s .

Remark 11 (benefit of the canonical representation): The canonical representation of (M1b) simplifies the shifting procedure since only the state vector is affected which, as a consequence, modifies the Hamiltonian and efforts variables. Any modification of the component law is similar, for analog electrical circuits, to swapping components on a PCB, as the connections remain identical, but the behaviour of the discrete components varies. Later on, the linearisation procedure also benefits from the canonical property of the interconnection.

We now establish a new model where the thermodynamical behaviour is linearised.

Model (M2): Linearised Thermodynamical Behaviour

Nonlinear acoustics relies mainly on two nonlinearities: the convection of matter, and the nonlinear state law relating the internal energy to the static pressure. Some physical systems, like the vocal tract, exhibit (small) acoustical fluctuations of the volumetric mass density, while their capacity to convect matter remains valid. In this case, it is safe to assume that the state law can be linearised while preserving the nonlinear convection law. This enables one to simplify the model and improve its numerical conditionning, as the compression part of the system becomes linear.

Here, we assume that the fluctuations of mass m in model (M1c) are small compared to its state at rest value m ‹ " m 0 such as

(App 1) : m m 0 " ρ ρ 0 ! 1 (5.16)
which is equivalent to assuming small fluctuations of ρ " ρ ´ρ0 , as we are considering a time-invariant volume V 0 (see (H0)).

To linearise the state law, we take the thermodynamical energy E c p mq of (M1c) in Eq. (5.14) and perform a series expansion of the second order so that

E c p mq " P 0 V 0 γ ´1 "ˆ1 `m m 0 ˙γ ´γ m m 0 ´1 " P 0 V 0 γ ´1 « p ¨¨¨1 `γ m m 0 `γpγ ´1q 2 ˆm m 0 ˙2q ¨¨¨¨γ m m 0 ´1ff `Op ˆm m 0 ˙3q « P 0 V 0 γ 2 ˆm m 0 ˙2 ,
where m is the (small) acoustical fluctuation of mass, such as m " m `Op m2 q.

Notation 4. Here, giving a generic physical variable a, â denotes its discretised counterpart (if needed). ã denotes its fluctuation around a state a rest value a 0 so that a " ã `a0 . a denotes the small fluctation so that a " a `a0 `Opã 2 q.

We define the state vector x 4 " r νL , νR , νT , ms where m is the acoustical fluctuation of mass around m 0 . The new Hamiltonian H 4 px 4 q is then quadratic w.r.t. m, such as

H 4 px 4 " r νL , νR , νT , msq " m `m0 2 2 0 p ν2 L `ν 2 R `ν L νR q `3pm `m0 q 2h 2 0 ν2 T `P0 V 0 γ 2 ˆm m 0 ˙2 .
The effort vector e 4 " ∇ x 4 H 4 px 4 q reads e 4 "

¨p m`m 0 q 2 0 p νL `ν R{2q p m`m 0 q 2 0 p νL{2 `ν R q 3p m`m 0 q νT{h 2 0 1 2 2 0 p ν2 L `ν 2 R `ν L νR q `3 2h 2 0 ν2 T `P0 V 0 γ m 2 0 m ‹ ‹ ‹ ‹ ‹ ‹ ' " ¨q L qR qT ă ψ ą Ω ‹ ‹ ‹ ‹ '
where only the term related to the state law P 0 V 0 γ m 2 0 m is linear and differs from the previous effort vector e 3 , showing that only the thermodynamical part of the model is affected by (App. 1).

Note that under (App. 1), the acoustical pressure inside the junction can now be retrieved by the well-known acoustical linear relation [CK13, Eq. (1.97)] m P 0 V 0{ρ 2 0 V 0 " c 2 0 ρ where c 0 "

a P 0 γ {ρ 0 .
Thanks to the canonical property obtained in (M1b) and inherited by (M1c), the differential algebraic remains similar, such as

¨dx 4 dt y 4 '" S 4 ¨e4 u 4 ',
where S 4 " S 3 , u 4 " u 3 and y 4 is the linearised version of y 3 under (App 1).

At last, we now establish a fully linear model where the phenomenon of mass convection is also neglected.

Model (M3): Linear Acoustics

A key hypothesis of linear acoustics is to neglect the impact of mass convection on the acoustical behaviour of the system [Fil+98, p. 12]. Some classical models of the vocal tract [Mae82; EL16] also make this hypothesis, as their physical setting is the one of linear acoustics.

Here, to neglect this phenomenon, we linearise the kinetic law of (M2) under the assumption that the fluctuations of the kinetic variables νi @i P tL, R, T u are small, which translates as

(App 2) :
νi " ν i `Op ν2 i q ! 1 @i P tL, R, T u.

(5.17)

This assumption impacts the kinetic part of the Hamiltonian, that we evaluate at m " m 0 ðñ m " 0, giving the following quadratic Hamiltonian H 5 px 5 q

H 5 px 5 " rν L , ν R , ν T , ms q " m 0 2 2 0 pν 2 L `ν2 R `νL ν R q `3m 0 2h 2 0 ν 2 T `P0 V 0 γ 2 ˆm m 0 ˙2
that we cast in the following matricial form

H 5 px 5 q " 1 2 x 5 » - - - - - m 0{ 2 0 m 0{2 2 0 0 0 m 0{2 2 0 m 0{ 2 0 0 0 0 0 3m 0{h 2 0 0 0 0 0 P 0 V 0 γ {m 2 0 fi ffi ffi ffi ffi fl Q x 5 implying e 5 " Qx 5
where Q is a matrix containing constant physical parameters.

The differential algebraic formulation is

p νL{2 `ν R q 3m 0 νT{h 2 0 ‹ ‹ ‹ ' (5.19)
where q is the acoustical mass rate and ψ the acoustical enthalpy.

Again, the acoustical pressure inside the cavity is given by mP 0 V 0{ρ 2 0 V 0 " c 2 0 ρ. At a given boundary S i i P tL, R, T u, the mass rate qi and enthalpy ψi are linked to the traditionnal acoustical pressure p i and flowrate U i by ˜pi

U i ¸" « ρ 0 0 0 1 {ρ 0 ff ˜ψ i qi ¸@i P tL, R, T u.
(5.20)

5.4 comparison of the junction models and state-space representation

In Tab. 2, we summarise the main differences and similarities for each junction model (M1a) to (M3). All the proposed models are passive and have a pH formulation.

Models (M1a-c) are equivalent: only their representation differ. (M1b-c) are reformulations of (M1a), allowing for: a pH canonical representation, in (M1b); and a nonlinear acoustical representation by defining the system as fluctuation around its state at rest that coincides with the atmospherical conditions, in (M1c). Models (M2) and (M3) incorporate approximations regarding the state law and the convection of mass. (M2) is equipped with a quadratic state law adapted for acoustical systems that exhibit mass fluctuations, but small variations of volumetric mass density and, therefore, pressure. (M3) is a fully linear acoustical system, where the state law and the mass convection phenomena are neglected. We analyse the state-space representation hereafter, to better understand the dynamics of our junction models. The linear state-space representation of model (M3) is obtained by setting x " x 5 , A " J 2 Q, B " G p2 , C " ´G p2 , D " 0 (see Eq. (5.13)), uptq " u 5 and yptq " y 5 (see Eq. (5.19)). The inputs and outputs are linked to traditional acoustical quantities by Eq. (5.20).

To retrieve the acoustical velocities and lumped volumetric mass density of the junction, one can invert the change of coordinate of Eq. (5.11), which gives M ´1xptq.

Remark 12 (Inverting the causality): As (M3) is linear, it is possible to invert the input/output causality so a given surface can be driven with an input mass rate (or flowrate) instead of an input enthalpy.

V 0 ρ{2 " v 2 mx `pv 2 c `v2 y q {3 ‰ m {2 r pν 2 L `ν2 R `νL ν R q { 2 0 `3ν T{h 2 0 s m`m 0{2 r p ν2 L `ν 2 R `ν L νR q { 2 0 `3 νT{h 2 0 s m`m 0{2 r p ν2 L `ν 2 R `ν L νR q { 2 0 `3 νT{h 2 0 s m 0{2 r p ν2 L `ν 2 R `ν L νR q { 2 0 `3 νT{h 2 0 s
H c pxq P 0 V 0{pγ´1q rp ρ {ρ 0 q γ ´γρ {ρ 0 s P 0 V 0{pγ´1q rp m {m 0 q γ ´γm {m 0 s P 0 V 0{pγ´1q rp m 0 `m{m 0 q γ ´γm {m 0 ´1s

Canonical J » - - - - - 0 0 0 ´3{ 0 V 0 0 0 0 ´3{h 0 V 0 0 0 0 0 3 { 0 V 0 3 {h 0 V 0 0 0 fi ffi ffi ffi ffi fl ´G p » - - ´1{2 0 ´3{2 0 0 0 1 {2 0 ´3{2 0 0 0 0 0 3 {h 0 0 fi ffi fl Model (M1) (M2) (M3) (H0-5)

Fluid Mechanics

Hpxq " H k pxq `Hc pxq analysis of eigenspaces We now investigate the dynamics of the model by looking at its eigenvectors. The characteristic polynomial of A is

P 0 V 0 γ {2p m{m 0 q 2 » - - 0 0 0 1 0 0 0 1 0 0 0 1 ´1 ´1 ´1 0 fi ffi fl » - 1 0 0 0 0 1 0 0 0 0 1 0 fi fl
λ 4 `P0 γ ρ 0 ˆ3 2 0 `h2 0 2 0 h 2 0 ˙λ2 " 0,
with a trivial root of multiplicity two for λ ‹ 1,2 " 0 and two complex conjugated solutions

λ ‹ 3,4 " ˘i c 0 ˜a3p 2 0 `h2 0 q 0 h 0 ¸,
where c 0 " b P 0 γ ρ 0 is the speed of sound. The associated eigenvectors are

v λ ‹ 1 " » - - - - - ´1 1 0 0 fi ffi ffi ffi ffi fl , v λ ‹ 2 " » - - - - - ´2 2 0{h0 0 1 0 fi ffi ffi ffi ffi fl , v λ ‹ 3 " » - - - - - Λ Λ Λ 1 fi ffi ffi ffi ffi fl and v λ ‹ 4 " » - - - - - ´Λ ´Λ ´Λ 1 fi ffi ffi ffi ffi fl where Λ " ip c 0{m 0 q a 2 0 h 2 0{3p 2 0 `h2 0 q.
The eigenvectors enables us to identify three main power exchange patterns. v λ ‹ 1 is a combination of νL and νR , which indicates a direct transfer of power from the left surface to the right surface, an axial incompressible flow. Likewise, v λ ‹ 2 is linked to an axial to transverse incompressible flow (left to the top surface). The last two remaining eigenvectors indicate a periodic exchange of power between the three fluid masses and the capacitor: the system acts as an oscillator where the capacitor acts as a spring and the fluid masses acts as the main inertia. This analysis is compatible with the equivalent electrical circuit depicted on Fig. 39.

As a corollary, we show how to compute the transfer matrix of (M3), a generalisation of the transfer function for multiple inputs multiple outputs (MIMO) systems.

Definition 3 (Transfer matrix [START_REF] Hinrichsen | Mathematical Systems Theory I : Modelling, State Space Analysis, Stability and Robustness[END_REF]): The transfer matrix Gpsq of a state-space system of form (5.21) is defined as Gpsq " D `Cps1 ´Aq ´1B.

(5.22)

For (M3), the transfer matrix is given by Gpsq " ´Gp5 ps1 ´J2 Qq ´1G p5 .

Remark 13 (straight bend geometry): Similary as in Ref. [START_REF] Chaigne | Acoustique des instruments de musique (2e {é}dition revue et augment{é}e)[END_REF]Sec. 7.7], a model of a straight bend (L shaped junction between two resonators) can be retrieved by setting q R " 0. On the equivalent electrical circuit, this is analogous to removing the component ν R .

summary and perspectives

In this chapter, we established five passive junction models and have shown how the use of a canonical formulation enables us to physically interpret the model (using equivalent electrical circuit), and simplifies the modelling process. We then examined the dynamics of the system by using a state-space representation, allowing for a better understanding of the power exchange patterns between the storing components.

The geometry we considered is simple, as well as the basis used for the projection of the continuous space-time equations. To account for more complex geometries and phenomena, our modelling methodology has to be generalised to account for a higher discretisation order. An interesting improvement would be to consider a geometry so that the side resonator is not orthogonal to the main bore. Finally, the implementation of this model in a pratical example would enable us to verify the intepretations we made. Fluid-structure interaction (FSI) problems arise when a deformable structure and a fluid are interacting together: the hemodynamic flow of blood in deformable vessels, the airflow at great speed around the flexible wing of an aircraft, interactions between a flag and a fluid, the glottal flow between the vocal folds. In these situations, the dynamics of the solid and the fluid phases are coupled, resulting in a multiphysical dynamical system. This coupling relies on three types of conditions [START_REF] Richter | Fluid-structure Interactions[END_REF]: the geometric conditions (i), the kinematic conditions (ii) and dynamical conditions (iii).

The geometric conditions (i) are a set of hypotheses that identify the different solid and fluid phases, and the volume they each occupy. Each phase is contained in a domain ω i ptq, possibly time varying, of boundary Bω i . The collection of regions (also called domains and subdomains) is a partition of the global time-invariant control volume Ω. At all times, the union of all regions ω i ptq perfectly covers1 the total region Ω, and the interior of the subdomains do not overlap.

As an example, we show on Fig. 40 a FSI dynamical system where a fluid phase (domain ω f ptq) is coupled with a solid phase (domain ω s ptqq. The coupling conditions at the common interior boundary Bω f X Bω s (see the red line on Fig. 40) are specified by the kinematic and dynamic conditions. ω s (t) At an interior interface and given an (H1) inviscid fluid assumption, the kinematic conditions (ii), i.e. the so-called normal velocity continuity, stipulate that the normal components of the fluid and solid velocity fields are equal, so that for all points in Bω f X Bω s , v f ¨n " v s ¨n, where v f is the velocity of the fluid phase, v s is the velocity of the solid phase and n is the outward-pointing normal vector to the considered interface.

ω f (t) ∂ω s ∩ ∂ω f v ∂
The dynamic conditions (iii) state the balance of forces at the boundary such as for all points in Bω f X Bω s , σ f n " σ s n, where σ f is the Cauchy stress tensor of the fluid phase, σ s is the stress tensor of the solid phase.

FSI problems are mostly tackled numerically. On the one hand, partitioned methods simulate each phase separately, using a Lagrangian approach for the solid, and a Eulerian approach for the fluid. The solid and fluid phases are then coupled using coupling algorithms2 . Partitioned methods are popular as they allow for the use of existing tools from structural and fluid mechanics. However, they do not always provide a satisfactory solution as they do not provide a unified description of the problem that encompasses the dynamics of each phase and the coupling conditions.

On the other hand, monolithic models aim at providing a unified description of the fully coupled system, the most successful one being the arbitrary lagragian eulerian (ALE) framework [START_REF] Donea | Arbitrary Lagrangian-Eulerian Methods[END_REF], also used to derive new partitioned methods.

Most of these methods do not guarantee the power balance of the model or the simulation, thus explaining the growing interest for pHs formulation of FSI problems, like the work of Cardoso Ribeiro in Refs. [Car16; CML19] and Mora Araque [START_REF] Alejandro | Port-Hamiltonian Modelling of Fluid-structure Interactions in a Longitudinal Domain[END_REF]. The pH framework is a well-suited formalism for the developement of energy-consistent models, because its modularity pairs well with the multiphysical nature of FSI problems. Some recent works aim at formulating the equations of fluid dynamics as a coordinate-free system [START_REF] Rashad | Port-hamiltonian Modeling of Ideal Fluid Flow, Part I: Foundations and Kinetic Energy[END_REF], in order to compute the interactions between the body of a flying bird and the surronding air [START_REF] Califano | Decoding and Realising Flapping Flight with Port-hamiltonian System Theory[END_REF]. However, these existing pH approaches are either not scalable [START_REF] Alejandro | Port-Hamiltonian Modelling of Fluid-structure Interactions in a Longitudinal Domain[END_REF] or do not account for time varying volumes [START_REF] Luiz | A Partitioned Finite Element Method for power-preserving discretization of open systems of conservation laws[END_REF].

To address this, Diagne and Maschke [START_REF] Diagne | Port Hamiltonian Formulation of a System of Two Conservation Laws with a Moving Interface[END_REF] provided an infinite-dimensional pH description of two one-dimensional systems coupled by a moving interface, where the time-varying domains occupied by each dynamical system are tracked by the use of indicator functions (also called color functions), a tool used in level-set methods [START_REF] Osher | The Level Set Methods and Dynamic Implicit Surfaces[END_REF] for instance. In computational fluid dynamics (CFD), the indicator functions are used to encode the geometric conditions [START_REF] Richter | Fluid-structure Interactions[END_REF]p. 262].

The modelling methodology proposed by Diagne and Maschke [START_REF] Diagne | Port Hamiltonian Formulation of a System of Two Conservation Laws with a Moving Interface[END_REF] 3 is threefold and described below:

Step 1 : the physical hypotheses and the geometric conditions (i) are stated;

Step 2 : the fully uncoupled pHs of the fluid and solid phases, and the associated time-varying subdomains is formulated; the domain ω s . Consider two color functions c f and c s , respectively representing ω f ptq and ω s ptq, and defined as

@i P s, f c i pr, tq " $ & % 1 if r P ω i ptq 0 if r R ω i ptq (6.1)
where r P R 3 is the position vector px, y, zq . Now, the domains ω f {s can be defined from the color functions as ω f ptq " tr P R 3 | c f pr, tq " 1u and ω s ptq " tr P R 3 | c s pr, tq " 1u. (6.2)

These color functions can vary with time and are associated with a PDE modelling the time-variation of the domain ω i ptq, so that c i is convected according to the following equation [START_REF] Mulder | Computing interface motion in compressible gas dynamics[END_REF] @i P ts, f u, Bc i Bt " ´gradpc i q ¨vB " v B ¨n, (6.3) where v B is the velocity of the boundary Bω f X Bω s and 0 otherwise, and where ´gradpcq generates the outward-pointing normal vector of the boundary n. When r P ω i , gradpc i q " 0, meaning that only the phenomena at the boundary have an effect on the dynamics of c i .

Remark 14 (inward-poiting normal vector): Due to the color function values, the gradient gradpc i q is postively oriented from the exterior of ω i (where c i " 0) towards the interior of ω i (where c i equals 1), and is orthogonal to Bω i .

To make sense of Eq. (6.3), we propose a new formulation based on distributions.

Reminder on Distribution Theory

We start by providing basic definitions and properties from distribution theory adapted from "Mathematics for the Physical Sciences" by Laurent Schwartz [START_REF] Schwartz | Méthodes mathématiques pour les sciences physiques[END_REF] (see also Ref. [START_REF] Schwartz | Théorie des distributions[END_REF]).

The common motionless domain Ω is a non-empty open set, equipped with a position vector r P Ω. We call the support of a real-valued function ϕ the smallest closed set, outside which ϕ " 0.

Definition 4 (Space of functions with bounded support, [START_REF] Schwartz | Méthodes mathématiques pour les sciences physiques[END_REF], Def. 1 p. 71):

The space DpΩq is the space of real-valued functions which are infinitely differentiable and have bounded supports on Ω:

DpΩq ": tϕ P C 8 c pΩqu . (6.4)

In the following, we sometimes abbreviate the notation DpΩq by D. Here, we consider test functions ϕ P DpΩq which vanish on the boundary of Ω. xT , ϕ 1 `ϕ2 y " xT , ϕ 1 y `xT , ϕ 2 y, (6.5a)

xT , λϕy " λxT , ϕy, @λ P R. (6.5b)

The distributions themselves form a vector space D 1 with the sum T 1 `T2 and the product λT defined by xT 1 `T2 , ϕy " xT 1 , ϕy `xT 2 , ϕy (6.6a)

xλT , ϕy " λxT , ϕy (6.6b) so that the scalar product xT , ϕy, for T P D 1 and ϕ P D so that the scalar product is a bilinear form.

The support of a distribution T [Sch65, p. 79] is the smallest closed set outside of which T is zero.

The following remark highlights a point that will be used later on in this chapter. We first consider the one-dimensional case where Ω :" ra, bs. We pick a function ϕ P DpΩq that vanishes at the boundary of BΩ. We use the integration by part procedure such as [Sch65, II, 2, 6, p. 80] A generalisation of this formula for a n-dimensional open set Ω Ă R n , n ě 2 is given below. Here, ϕ is a test function of vectorial values on Ω. The derivative of T , denoted by its gradient gradpT q, reads x gradpT q , ϕy " Then, the first order derivative of T is given by the following property.

Property 2 (Gradient of a distribution): The gradient of a distribution T P D 1 pΩq is defined as @ϕ P DpΩq, x gradpT q , ϕy " ´xT , divpϕqy , (6.7)

where ϕ is a vectorial test function in DpΩq. Similarily, δ r evaluates ϕ at the point r P Ω, that is, xδ r , ϕy " ϕprq.

A
We now proceed to the definition of the color distribution, a term that replaces the term color function to avoid any ambiguity.

Color Distributions

We start by specifying the geometric conditions for the general case of the full FSI problem. Consider the domain Ω (global time-invariant control volume), a nonempty connected open set, split into N ω time varying open subdomains ω i ptq associated with each fluid and solid phase, as depicted on Fig. 41.

The collection ω iP 1, Nω ptq defines a cover of Ω [GG99, p. 19] the union of which is the whole region Ω, such as

Ω " Nω ď i"1 ω i ptq.
(6.8)

The intersection of two contiguous regions is the common interface ω i X ω j " Bω i X Bω j @ ti, j P 1, N ω | i ‰ ju , (6.9) stating that the fluid and solid phases never overlap, but share common interfaces.

If ω i and ω j are not contiguous, their intersection is the empty set. Remark 16 (equality of Ω and the union of its cover): The sets Ω and Ť Nω i"1 ω i ptq are equal in the sense of measure theory, meaning that Ωz

Ť Nω i"1 ω i ptq is of null measure.
The global domain Ω is equipped with the scalar product x¨, ¨y. Moreover, for each member of the cover and each t, ω i is equipped with its own scalar product x¨, ¨yω i ptq restricted to the domain ω i ptq defined as xa , by ω i ptq " ţ ω i ptq ab dV . We denote by DpΩq the set of infinitely differentiable functions with a compact support on Ω and by D 1 pΩq the vector space of distributions on Ω.

We can now define the concept of color distribution.

Definition 6 (color distribution): Let t Þ Ñ ω i ptq be a continuously time-varying region in Ω. Then, we define the color distribution of ω i as c i p¨, tq P D 1 pΩq such that @ t, @ ϕ P DpΩq, xc i p¨, tq , ϕy " x1 , ϕy ω i ptq , (6.10)

where 1 is the identity distribution on Ω.

Thus, a color distribution assigns a function ϕ to the volume-integral of ϕ on the domain ω i ptq. In the following, we omit the variable t.

Property 3 (summing color distributions): Let tc i u, i P 1, N ω be a (static) cover of Ω. The sum of all the color distributions c i applied to ϕ P DpΩq yields

Nω ÿ i"1 xc i , ϕy " Nω ÿ i"1 x1 , ϕy ω i " x1 , ϕy.
(6.11)

A corollary is that this property is also true for any time-varying fonction cp¨, tq P DpΩq, in the sense that, for all t,

Nω ÿ i"1 xc i p¨, tq , ϕy " Nω ÿ i"1
x1 , ϕy ω i ptq " x1 , ϕy, (6.12) because the difference between Ω and Ť Nω i"1 ω i ptq is of null measure (see Rk. 16).

gradient of a color distribution -We can now define the spatial derivative of c i (see Eq. (6.3)), that is, the distribution gradpc i p¨, tqq, denoted by gradpc i q in the following. More generally and for the sake of conciseness, we omit the space-time dependence in the following when it is not ambiguous. As a consequence of Property 2, applying gradpc i q to a function ϕ P DpΩq computes the flux of ϕ through the boundary of ω i . Using property 2 yields

x gradpc i q , ϕy (6.7)

" ´xc i , divpϕqy, which is equivalent, under Def. 6, to the integral on the subdomain ω i ptq Ă Ω ´xc i , divpϕqy " ´x1 , divpϕqy ω i .

Finally, using the divergence theorem [LSS09, Th. 6.1] yields the result ´x1 , divpϕqy ω i " ´¡ ω i divpϕq dV " ´ij Bω i ϕ ¨n dS.

We may now state the full definition.

Definition 7 (Gradient of a color distribution): Let ω i be a domain in Ω, associated with a color distribution c i P D 1 pΩq. Then, the distribution gradpc i q assigns a function ϕ P DpΩq to the incoming flux of ϕ through the boundary of ω i , denotes by Bω i , so that

x gradpc i q , ϕy " ´ij Bω i ϕ ¨n dS, (6.13)

where n is the outward-pointing normal vector of the boundary Bω i .

This definition is illustrated on Fig. 42 for a generic subdomain ω i . The gradient of a distribution gradpc i q can be seen as the surface generalisation of the one dimensional Dirac (or delta) function at the boundary of the associated domain. 6.3), within the distribution framework. To this end, we now work with test function in DpΩ ˆT q. We define by T the arbitrary compact time interval rt 1 , t 2 s with t 1 ă t 2 . We pick a compact and infinitely differentiable test function ϕpr, tq in the space-time domain DpΩ ˆT q, which is zero on BpΩ ˆT q. where the first term vanishes as the ϕ is zero on BT . In the above equation, v B is the velocity of the boundaries, defined on Ş Nω i"1 ω i . By definition, v B is null on the boundary of the (motionless) domain Ω.

We transform the surface integral in a volume integral on Ω by applying the divergence theorem [START_REF] Yamaguchi | Engineering Fluid Mechanics[END_REF] (in a general form that accounts for a scalar, vectorial, or tensorial quantity), so that (omittin the space-time variables)

ż T ij Bω ϕ pv B ¨nq dS dt " ż T x1 , div `ϕv B ˘yω i dt " ż T xc i , div `ϕv B ˘y dt " ż T ¡ Ω div `ci ϕv B ˘´`g radpc i q ¨vB ˘ϕ dV dt " ż T ¨ij BΩ `ci ϕv B ˘¨n dS 0 ‹ ‹ ‹ ‹ ' dt ´żT ¨¡ Ω `gradpc i q ¨vB ˘ϕ dV ' dt
where the first term is zero because v B is null on BΩ (boundary of the motionless global control volume).

Finally, we obtain the following definition ż T ˆx Bc i Bt , ϕy `x gradpc i q ¨vB , ϕy ˙dt " 0 (6.14) which states, under Def. 7 and in the sense of distributions, that the time variation of the domain ω i is equal to the normal component of the velocity of its boundary. This is equivalent to Eq. (3.3) in Ref. [START_REF] Diagne | Port Hamiltonian Formulation of a System of Two Conservation Laws with a Moving Interface[END_REF].

Definition 8 (Time derivative of a color distribution): Let t Þ Ñ ω i ptq be a continuously time-varying region in Ω of boundary velocity v B . Let c i P D 1 pΩ ˆT q be the color distribution associated with ω i . Then, the distribution Bc i Bt assigns a function ϕ P DpΩ ˆT q to the flux induced by the motion of the subdomain ω i at its boundary Bω i , such as

x Bc i Bt , ϕy ΩˆT " ´x gradpc i q ¨vB , ϕy ΩˆT " ż T ¨ij Bω i ϕ pv B ¨nq dS ' dt (6.15)
where x ' , ' y ΩˆT is the natural scalair product on tΩ ˆT u, and where n is the outward-pointing normal vector to the boundary Bω i .

Extrapolation of Physical Variables on the Common Domain

Consider a generic physical quantity p i pr, tq that is C 1 regular on the subdomain ω i (p i P C 1 pA i , R n q with A i " tω i ptq ˆT , @tu Ď Ω ˆT ) which can be scalar (e.g. internal energy, volumetric mass density), vectorial (e.g. velocity) or tensorial (e.g. stress and strain tensors). We denote pi P C 1 P W pΩ ˆT q (the space of piecewise smooth function on Ω ˆT ) the extrapolation of p i over Ω ˆT ( pi is equal to p i on A i and zero elsewhere) so that xc i pi , ϕy ΩˆT " xc i p , ϕy ΩˆT (6.16)

where we introduced the extrapolated variable p P C 1 P W pΩ ˆT , R n q. Here, pi takes the value p i on ω i , and 0 outside ω i (that is pi " c i pi ).

We generalise this definition to account for a collection of physical variables pi of the same nature (e.g. mass density, velocity). Consider the sum of each physical variables p i of the N ω subdomains tω i u iP 1, Nω so that Nω ÿ i"1 xc i pi , ϕy ΩˆT " x Nω ÿ i"1 c i pi , ϕy ΩˆT " x p , ϕy ΩˆT using the bilinearity property of the scalar product and remark 15.

Here, the extrapolated variable p P C 1 P W pΩ ˆT , R n q takes the value p i when r P ω i ptq @i P 1, N ω . Definition 9 (extrapolated variable): Consider a collection of N ω physical variables pi P P, with P " C 1 P W pΩ ˆT , R n q, of the same nature. We define the extrapolated variable ppr, tq P P that takes the value pi when r P ω i ptq such as @ϕ P DpΩ ˆT q, x Nω ÿ i"1 c i pi , ϕy ΩˆT " x p , ϕy ΩˆT , i.e. formally p "

Nω ÿ i"1 c i pi .
Property 4 (Applying a color distribution on an extrapolated variable): Applying the color distribution c i , associated with the time-varying subdomain ω i ptq, to an extrapolated variable p P P yields the variable pi , that is @ϕ P DpΩq, xc i p , ϕy ΩˆT " x pi , ϕy ΩˆT Corollary 1 (Empty intersection of domains): Under Eq. (6.9) and remark 15, applying a color distribution c i P D 1 pΩ ˆT q to a variable pj P P yields 0 for i ‰ j, that is, @i, j P 1, N ω , @ϕ P DpΩ ˆT q, i ‰ j, ùñ xc i p j , ϕy ΩˆT " 0 , because the support of c i and p j do not overlap (see Rk. 16).

color-distribution-augmented multiphase problem

In this section, we introduce the fluid pH model that we consider. Then, we detail the methodology to derive pHs augmented with color distributions. To generalise the procedure for a generic infinite-dimensional pHs, see the approach of Ref. [START_REF] Diagne | Port Hamiltonian Formulation of a System of Two Conservation Laws with a Moving Interface[END_REF].

We consider N f fluid phases. Each fluid phase is indexed by i P F, F being a n-tuplet containing the index of each fluid phase. We first isolate a single fluid phase contained in the subdomain ω i @i P F. fluid model -We assume that the i-th fluid phase verifies the hypotheses (H1, 3-5). The fluid and flow are described by the velocity v i pr, tq, the volumetric mass ρ i pr, tq for r P ω i , and the system of conservation equations [START_REF] Leal | Advanced Transport Phenomena: Fluid Mechanics and Convective Transport Processes[END_REF] describing the time variation of mass Bρ i Bt `divpρ i v i q " 0 (6.17) and the time variation of momentum

Bρ i v i Bt `divpρ i v i b v i q `gradpp i q " 0. (6.18)
The energy stored by the fluid is the sum of the kinetic and internal energy

E f " ¡ ω ˆ1 2 ρ i v i ¨vi `ρi i pρ i q ˙dV (6.19)
where i pρ i q is the specific internal energy, related to the pressure by [vM02]

p i " ρ 2 i B pρ i q Bρ i . (6.20)
port-hamiltonian form -We consider the model of Ref. [START_REF] Alejandro | Port-Hamiltonian Modelling of Fluid-structure Interactions in a Longitudinal Domain[END_REF] under hypotheses (H1, 3-5). Choosing v i and ρ i as energy variables gives the following Hamiltonian

H f ,i pX f ,i " " v i , ρ i ‰ q " ¡ ω i H f ,i pX f ,i q dV with H f ,i pX f ,i q " 1 2 ρ i v i ¨vi `ρi pρ i q leading to the efforts e i " δH f ,i pX f ,i q δX f ,i " « δHpv i q δv i δHpρ i q δρ i ff " « ρ i v i 1 2 v i ¨vi ` pρ i q `ρi B pρ i q Bρ ff " « π i ψ i ff (6.21)
where π i is the volumetric momentum and ψ i the total specific enthalpy.

The dynamics of the pHs is given by

B Bt ˜vi ρ i ¸" « ´Gω,i{ρ i ´gradp ¨q ´divp ¨q 0 ff « π i ψ i ff (6.22)
where divpq and gradpq are formally skew-adjoint according to Ref. [START_REF] Zwart | The Linear Wave Equation on N-dimensional Spatial Domains[END_REF] and where G ω,i , first introduced in Ref [Car16, 6.2, p.

-73], is a skew-symmetric matrix called a gyroscopic operator, that is responsible of the power transfers between the different components of the velocity due to the vorticity of the flow. It is defined as

G ω,i " ¨0 ´ωz,i ω y,i ω z,i 0 ´ωx,i ´ωy,i ω x,i 0 ‹ ‹ ' with ω i " rot v i " » - - - ω x,i ω y,i ω z,i fi ffi ffi fl .
(6.23)

The power-balance reads dH i pX i q dt " ´ij Bω i pπ i ψ i q ¨n dS. (6.24)

augmented state vector and hamiltonian -Following the procedure of Ref. [START_REF] Diagne | Port Hamiltonian Formulation of a System of Two Conservation Laws with a Moving Interface[END_REF], we define the state vector of all the fluid pHs comprised of the extrapolated variables ṽ and ρ, and the associated color distributions such as X f " r ṽ , ρ, c F p1q , . . . , c F pN f q s . The total Hamiltonian is the sum of each energy density H i accross each fluid phase i P F such as

H f pX f q " ÿ iPF H i pX q " ÿ iPF ¡ ω i H f ,i pX f q dV (6.25)
The effort variables are given by the variational derivative of the Hamiltonian w.r.t. the state vector. We first take the component associated with ṽ as an example, so that

e ṽ " δH f pX f q δ ṽ " δ δ ṽ ¨ÿ iPF ¡ Ω c i H f , i pX f q dV ' " ÿ iPF δ δv i ´ÿ iPF ¡ ω i H f , i pX f q dV ¯Bv i B ṽ c i " ÿ iPF ˆδH f ,i δv i (6.21) " π i c i " π
using the chain rule for the variational derivative [GR96, p.-37, Eq. 8]. Likewise, we have e ρ " ψ.

The effort e c j , associated with the j-th color distribution is given by

@j P F, e c j " δH f ,i pX f q δc j " δ δc j ¨ÿ iPF ¡ Ω c i H f ,i pX f q dV '" H f , j pX f q.
Then, the effort vector e f of the fluid phases, gathering every effort variable, reads

e f " δH f ,i pX f ,i q δX f ,i " " π , ψ, H f , F p1q pX f q, . . . , H f F pN f q pX f q ı . (6.26)
We now establish the conservation equation for each component of the state vector, starting with the volumetric mass density. conservation of mass -For the dynamics of ρ, we give the following theorem.

Theorem 1 (dynamics of ρ in the sense of distribution): In the sense of distributions in Ω, we have where we used the linearity property of the distribution. We then write d dt x ρi c i , ϕy in an integral form (omitting variable t)

ÿ iPF d dt x ρi c i , ϕy " ÿ iPF d dt ¡ Ω c i ρi ϕ dV " ÿ iPF d dt ¡ ω i ptq ρi ϕ dV .
Then, applying the Leibniz-Reynolds theorem, it follows that

ÿ iPF d dt ¡ ω i ptq ρi ϕ dV " ÿ iPF ¡ ω i ptq B ρi Bt ϕ dV 1 `ÿ iPF ij Bω i ptq ρ i v B ¨n ϕ dS. 2
The partial derivative of ρi w.r. We write the dynamics of ρ in a local form (omitting the space-time variables)

B ρ Bt `divp πq " ÿ iPF DIV c i π 1 ´ÿ iPF Γ c i π 2 (6.32)
where we highlighted: in blue , on the left-hand side, the terms accounting for the dynamics within each domain ω i , where the distribution DIV c i and Γ c i are null; in orange , on the right-hand side, the terms related to the boundary interactions, on Bω i where the distribution DIV c i and Γ c i restore the jump conditions with the neighbouring phases.

We now examine at the right-hand side of Eq. (6.32). Under theorem 1, it defines the effective mass-flow at a discontinuity [Fil+98, Eq. (1.18)] (i.e. the boundaries Bω i ):

1 ´2 " ´ÿ iPF ij Bω i pπ i ¨nq dS `ÿ iPF ij Bω i pρ i v B ¨nq dS " ÿ iPF ij Bω i pρ i v B ´πi q ¨n dS.
Both terms combine depending on the nature of the boundary. In the particular case of an impermeable boundary, the kinematic condition gives v i ¨n " v B ¨n, and the mass flow is null. As of now, we did not make any assumptions regarding the nature of the boundaries Bω i . conservation of momentum -Applying theorem 1 on d dt x ṽ , ϕy, ṽ P P and ϕ P DpΩq yields a similar result,

x B ṽ Bt `grad `ψ ˘`G ω π ρ , ϕy " ÿ iPF DIV c i ψ 1 ´ÿ iPF Γ c i ψ 2 (6.33)
where term 1 is the input enthalpy at the boundary, term 2 is the enthalpy induced by the displacement of the interface.

color distribution -The dynamics of the i-th color distribution is given by Eq. (6.14)

@i P F, Bc i Bt " ´vB ¨gradpc i q 2 . (6.34)
We gather Eqs. (6.32) (6.33) and (6.34) below

B ṽ Bt " ´grad `ψ ˘´G ω π ρ `ÿ iPF DIV c i ψ ´ÿ iPF Γ c i ψ (6.35a) B ρ Bt " ´divp πq `ÿ iPF DIV c i π ´ÿ iPF Γ c i π (6.35b) @i P F, Bc i Bt " ´gradpc i q ¨vB . (6.35c)
which governes the dynamics of the uncoupled multiphase flow dynamical system. Note that with the port-Hamiltonian notations we provided, it is possible to identify the components of the effort vector in the dynamics of ρ and ṽ. In Eqs. (6.35), we identify connecting ports indicated by the presence of the boundary distributions DIV c i and Γ c i . The velocity of the boundary v B acts on the dynamics of the phases through the distribution Γ c i . Depending on the physics at hand, this velocity may be driven by a neighbouring fluid phase, or by a solid phase (the formulation of which has not been addressed in this thesis).

summary and perspectives: towards a port-hamiltonian formulation

In this chapter, we gave a new definition of color functions within the scope of distributions. To this end, we developped tools and definitions to finally express the dynamics of multiple phases of fluid in the sense of distributions. We defined boundary distributions, that map the interaction phenomena at the boundaries to the dynamics inside the phases (contained in ω i ). We introduced the port-Hamiltonian formulation of the fluid phases, in order to inspire forthcoming studies to propose a full port-Hamiltonian formulation augmented with color distributions.

towards a port-hamiltonian formulation -Inspired by the original work of Diagne and Maschke [START_REF] Diagne | Port Hamiltonian Formulation of a System of Two Conservation Laws with a Moving Interface[END_REF], we give some elements that could lead to the formal definition of port-Hamiltonian systems that include color distributions.

Consider the N f fluid phases, as defined in the previous section. We propose a possible expression of the dynamics of the fluid phases that includes: a linear Hamiltonian operator J ; abstract "boundary operators" B f (gathering terms associated with DIV c i in Eqs. (6.35)) and D f (gathering the terms associated to Γ c i ).

On the one hand, B f maps the interactions with the neighbouring domains (jumps) to the internal dynamics of the considered phase. On the other hand, D f maps the input velocity of the boundary v B to the effect induced by the displacement of the boundaries on the dynamic of the system. Thus, it could be possible to express, as in Ref. [START_REF] Diagne | Port Hamiltonian Formulation of a System of Two Conservation Laws with a Moving Interface[END_REF], the dynamics of the multiphase flow in the following manner: We show a possible compact formulation that includes: the internal dynamics within each subdomain, highlighted in blue ; the terms related to the interactions at the boundaries mapped to the internal dynamics by the means of B f and D f , highlighted in orange ; and the power-conjugated ouputs f f B and p B . It is similar to the algebraic-differential formulation for the finite-dimensional pHs, and reads ¨Bt ṽ

BX f Bt " J f e f `Bf pX f qe B f ´Df pX f qv B where J f " » - - - ´G
B t ρ B t c F p1q B t c F p2q . . . B t c F pN f q f f B p B ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' " » - - - - - - - - - - - - - - - - - J f B f D f O pN f ˆ1q B ‹ f O p1ˆN f q O p3ˆ3q D ‹ f fi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi fl ¨δH δ ṽ " π δH δ ρ " ψ δH δc F p1q " H 1 δH δc F p2q " H 2 . . . δH δc F pN f q " H N d e f B v B ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ '
where B ‹ f and D ‹ f denote the dual of B f and D f .

perspectives -Establishing the color-distribution-augmented system of the solid phases will require the use of a more complex mathematical setting. Indeed, solid mechanics rely on vectorial and tensorial quantities. The investigation of a possible port-Hamiltonian structure, as in Ref. [START_REF] Diagne | Port Hamiltonian Formulation of a System of Two Conservation Laws with a Moving Interface[END_REF], should lead to a well posed and physically interpretable formulation of FSI problems.

G E N E R A L C O N C L U S I O N

contributions

In this thesis, we proposed modelling tools and passive physical models of the vocal apparatus, with a particular focus on the vocal tract. The contributions of the thesis are listed below, sorted by chapter.

port-hamiltonian systems -After introducing this modelling framework, we proposed, in Sec. 2.3.3 an improved version of the discrete gradient method less prone to cancelation errors. Then, we introduced a new procedure to formulate a shifted pHs, the physical quantities of which are written as fluctuations around a chosen state at rest.

power-balanced fsi lumped-parameter modelling of the vocal tract -A new power-balanced FSI model of the vocal tract has been established. Starting from a set of geometrical and physical hypotheses, we proposed a pH formulation of a model of a single tract. We interpreted and compared it to an existing model of the literature using equivalent electrical circuits. To build the full vocal tract, we proposed, in Sec. 3.2, a general constraint-based assembly method that accounts for an arbitrary configuration setting (e.g. simple connection, branching of multiple resonators). Then, in Sec. 3.2.4, we applied this method to the vocal tract, and showed how taking advantage of the structure of a pHs simplifies the resulting constrained model.

numerical experiments -We performed simple numerical experiments that confirmed the expected capacities of the tract model, namely: its capacity to convect matter, to account for the movement of the walls, to exhibit acoustical resonances, and to be passive. Then, using a very simple glottal source, we simulated a simple coarticulation.

a passive three-port junction for fluid mechanics and acoustics -We illustrated our modelling methodology by establishing a passive threeport junction model, form which we derived five other models. Each model is equipped with particular hypotheses that are targeted at: modelling fluid mechanics, nonlinear acoustics, semi-nonlinear acoustics and linear acoustics. Tab. 2 recapitulates the similarities and differences between each junction model. The use of a canonical interconnection allowed for the physical interpretation of the nonlinear acoustical model. We linked our work to classical linear acoustical approaches using the state space representation of the fully linear junction.

an infinite dimensional ph formulation for fsi problems using color distributions -This chapter is a preliminary work on the use of indicator functions (also called color functions) within the infinite-dimensional pH framework. We proposed a new definition of color functions that relies on distributions. We built an uncoupled formulation of a generic multiphase fluid problem. This formulation accounts for the dynamics of each fluid phase, within their respective domains, and is well suited to specify the interactions at the boundary between each phase.

perspectives

port-hamiltonian systems -The improvements made on the discrete gradient methods were carried out without looking at the overall implicit time-discrete problem. Taking the rate of convergence of the Newton-Raphson solver into account could lead to a smarter form of the discrete gradient. Taking constraints into account would allow for simulation of DA-pHs. At last, generalising the quadratisation method [START_REF] Lopes | Approche passive pour la modélisation, la simulation et l'étude d'un banc de test robotisé pour les instruments de type cuivre[END_REF] to account for constraints could lead to simulations compatible with real time. Moreover, using the power-balanced time finite element methods of arbitrary accuracy order developed by Müller [START_REF] Müller | Time-continuous Power-balanced Simulation of Nonlinear Audio Circuits : Realtime Processing Framework and Aliasing Rejection[END_REF] could be an interesting perspective, as it would enable us to reduce the sampling frequency since the interpolating functions are not band-limited.

Regarding the shift method, we only considered pHs equipped with a canonical interconnexion matrix. An important work is to generalise the method for all pHs and study its relevance on the numerical conditionning of the discrete-time problem. Generalising the method to shift a system around a trajectory would allow for interesting developpements regarding the study of bifurcations.

power-balanced fsi lumped-parameter modelling of the vocal tract -The proposed vocal tract model will be a useful tool to investigate the interpretation strategies of wind instruments players. Taking the energetic contributions of the vocal tract into account will help to understand its role in this situation. The geometry of our tract model is cartesian, and only the cross section varies with time. Considering cylindrical or elliptical geometries with a variable length would yield a better modelling of the fluid mechanical phenomena and allow for more realistic synthesis of the vowels, the associated geometry of which requires the vocal tract to lengthen (e.g. vowel [u]). Adding turbulence losses and noise generation will also help and improve both the realism and numerical conditioning of the tract model. numerical experiments -Using a numerical method tailored for constrained systems would allow for more stable and realistic simulations. Adding a switchable pHs model of a tract will allow for the complete closure of the vocal tract and the synthesis of stop consonants. Coupling the vocal tract with an articulatory model will allow for more complex and realistic coarticulation settings. a passive three-port junction for fluid mechanics and acoustics -Using this junction in a real setting (e.g. to model the lateral holes on wind instruments) would give more information about the realism and relevance of this passive model. In the thesis, we provided some elements to allow for a higherorder discretisation, which could lead to accounting for more complex geometries (e.g. cylindrical or conical geometries) and high frequency phenomena.

an infinite dimensional ph formulation for fsi problems using color distributions -This first work should serve as an inspiration to build a thorough theory, relevant for FSI problems. We considered the particular physical model of a perfect fluid. Generalising this method to any infinite-dimensional pHs and including dissipation phenomena could be useful for other applications. As an example, consider the interaction of a bow with a string, where the contact surface (i.e. where the hair and string meet) varies with time. The coupling and time-tracking of the contact surface may be naturally described (in a modular way) using color distributions. example: the mka system We choose X s " diag pπ s , x ks q, where we choose π s " 4.1 kg ¨m ¨s´1 and x ks " 0.1 m as the maximum value of both quantities on the numerical experiment shown on Fig. 7. We apply the adimensionalisation method and perform the same exact numerical simulations which results are plotted on Fig. 44. The total internal energy E comp is derived using the hypothesis (H3-5). The change of internal energy pT q, function of the temperature T , is governed by Joule's law, the perfect gas law 1 and the adiabatic process equation d " C V dT ; T " p {ρr 0 and p " P 0 p ρ {ρ 0 q γ (B.1.1)

where C V is the heat capacity at constant volume, r 0 the specific universal gas constant, and ρ 0 the volumetric density at rest. Integrating Joule's law between T 0 and T yields pT q ´ pT 0 q " C V pT ´T0 q. We substitute the perfect gas law and the adiabatic process equation in place of T and T 0 , assuming a rest state T 0 , ρ 0 and P 0 . Finally, we choose to set the remaining constant pT 0 q at ´P0{ρ 0 2 . The volumetric compression energy density U pρq " ρ pρq then reads

U pρq " P 0 γ ´1 "ˆρ ρ 0 ˙γ ´γ ρ ρ 0  (B.1.2)
where γ " r 0{C V `1 is the heat capacity ratio. The last step is to take (H7) into account (small fluctuation of ρ around its rest state), which mathematically translates to ρ´ρ 0{ρ 0 " ρ{ρ 0 ăă 1. Using this, we perform a second order series expansion 1 Shown here in its specific form. 2 This constant is usually arbitrarily set to zero. Here, we wish for the energy to be compatible with linear acoustics. It sets the minimum of the energy at the rest state T 0 , ρ 0 , P 0 rather at the absolute zero, which is consistent with the assumptions of linear acoustics (H7), as the heat capacific at constant volume of air is more or less constant around atmospheric conditions (H5).

i modelling approach and tools 1

 1 Figure 1Front cut of the vocal apparatus. . . . . . . . . . . . . . . . . . . . . . . Figure2Front cut of the larynx . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure3A vocal tract approximated by concatenated cylinders. . . . . . . . . . . Figure4A linear mass-spring-damper (MKA) system. . . . . . . . . . . . . . . . Figure5A two-mass-spring-damper system where m 1 and m 2 are rigidly connected. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure6A one-dimensional acoustical pipe with rigid and motionless walls. . . . Figure7Numerical simulation results for the mass-spring-damper system. . . . . Figure8Numerical results for different version of the discrete gradient (DG) method. Figure9Example of the shift procedure for a one-dimensional Hamiltonian function. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 10 Simplified sagittal cut of the vocal tract . . . . . . . . . . . . . . . . . . Figure 11 The lumped-parameter approach. . . . . . . . . . . . . . . . . . . . . . . Figure 12 Considered junction between K tracts through a common interface . . . Figure 13 2D cut of a tract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 14 Velocity field decomposition. . . . . . . . . . . . . . . . . . . . . . . . . . Figure 15 Equivalent electrical circuit of tract model . . . . . . . . . . . . . . . . . Figure 16Equivalent electrical circuit of the tract model of Maeda . . . . . . . . . Figure17A full vocal tract shown as a network of interconnected tracts . . . . . . Figure18A basic two-tract example. . . . . . . . . . . . . . . . . . . . . . . . . . . Figure19Oriented graph of the network in Fig.17. . . . . . . . . . . . . . . . . . Figure 20 Example: directed graph of the two-tract case . . . . . . . . . . . . . . . Figure 21 Example of a vocal tract with N t " 5. . . . . . . . . . . . . . . . . . . . . Figure 22 Directed graph of example 2 . . . . . . . . . . . . . . . . . . . . . . . . . Figure 23 Spring-Damper system . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 24 Equivalent electrical circuit of the shifted tract model augmented with simple resistors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure 25 Two-tract system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Figure26Height and total mass versus time at the dynamical phase. . . . . . . . . Figure27Local error on the powerbalance for the two-tract model. . . . . . . . . . Figure28Frequency responses of the two-tract model for both static phases . . . . Figure29Frequency response of a uniform, closed-open, 17-centimeter-long tube discretised with N t " 1, 22 tracts. . . . . . . . . . . . . . . . . . . . . Figure30Simulated geometry of the vocal tract. . . . . . . . . . . . . . . . . . . . Figure31Frequency response of the vocal tract for phases (1) and (3). . . . . . . Figure32Power flows of each tract during step (1) and (2). . . . . . . . . . . . . . Figure33Flow rate at the lips and its time derivative for a simple coarticulation. . Figure34Numerical results for the coarticulation a-to-i . . . . . . . . . . . . . . . Figure35Input enthalpy and output mass rate at both extermities of the vocal tract. Figure36A generic three port junction between 3 systems. . . . . . . . . . . . . . Figure37The five different junction models . . . . . . . . . . . . . . . . . . . . . .

  example for the MKA using PyPHS . . . . . . . . . . . . 119 xix A C R O N Y M S ALE arbitrary lagragian eulerian CFD computational fluid dynamics DA-pHs differential-algebraic port-Hamiltonian system DA-pH Differential-algebraic port-Hamiltonian DAE differential-algebraic equation DG discrete gradient DoF degrees of freedom FEM finite element method FSI fluid-structure interaction IDG improved discrete gradient MIMO multiple inputs multiple outputs ODE ordinary differential equation PCB printed circuit board PDE partial differential equation PFEM partitioned finite element method pHs port-Hamiltonian sytems pH port-Hamiltonian xx I N T R O D U C T I O N context

Figure 1 :

 1 Figure 1: Front cut of the vocal apparatus. In red, the respiratory group. In blue, the phonatory group. In orange, the articulatory group. From the illustrations of Anatomie pour la voix by Blandine CALAIS-GERMAIN © 2013 -ADVERBUM for the editor D ÉSIRIS, GAP -FRANCE.

Figure 2 :

 2 Figure 2: Front cut of the larynx, from the illustrations of Respiration by Blandine CALAIS-GERMAIN © 2005 -ADVERBUM for the editor D ÉSIRIS, GAP -FRANCE.

  Figure 4: A linear mass-springdamper (MKA) system.

velocity of the mass m 1

 1 reaction force of the spring reaction force of the dashpot force received by the exterior (2.6)

Figure 5 :

 5 Figure5: A two-mass-spring-damper system where m 1 and m 2 are rigidly connected.

Figure 7 :

 7 Figure7: Simple simulation of a MKA system performed with the PyPHS library. Here: m " 0.5 kg, k " 4441 N ¨m´1 so that f " 15 Hz; a " 5 kg ¨s´1 . On the top, the time evolution of the energy variables. In the middle the time evolution of the input F ext . On the bottom left, the phase space, and on the bottom right the power balance.

  Eq. (1.11) and Eq. (1.12) with DG = √ DG, Eq. (1.11) IDG, Eq. (1.12) (b) IDG versus the classical DG method.

Figure 8 :

 8 Figure 8: Numerical results for different version of the DG method. Here the Hamiltonian is k sinhpxq with k " 1.

Figure 10 :

 10 Figure 10: Simplified sagittal cut of the vocal tract with its main cavities (blue: pharynx, green: nasal cavity, orange: oral cavity); the nasal branching occurs where the three cavities meet.

Fluid

  Continuous 2D FSI problem and space discretisation into tracts.

Figure 11 :

 11 Figure 11: Lumped-Parameter approach. At each interface, the connection variables are: the mass rate q and enthalpy ψ for the open surfaces; velocity v and force F for the impervious surface (see Sec. 3.1.1 below).

Figure 13 :

 13 Figure 13: 2D cut of a tract; in red, the port variables; the grey arrows denote the main direction of the flow; S L , S R are open and motionless surfaces; S B is impervious and motionless; S W is impermeable and move with uniform velocity v W .

  accounts respectively for the uniform axial flow, the incompressible pumping of the moving wall and the change of velocity due to the compressibility of the fluid, as illustrated on Fig.14. The boundary condition at the moving wall (3.2) now reads v y ptq " v W ptq.

Figure 14 :

 14 Figure 14: Velocity field decomposition, from left to right: uniform axial flow; pumping effect due to the moving wall; compressibility effect. Here, we assume that v mx , v y and v c are positive.

Figure 15 :

 15 Figure 15: Equivalent electrical circuit of tract model with the storing components (blue), the memoryless components (green) and the external ports (red) displayed with the same topology as on Fig. 13 (left, right and top surfaces); on the left, analogy used for each part of the system (see Tab. 1).

Figure 16 :

 16 Figure 16: Equivalent electrical circuit of the tract model of Maeda [Mae82] under the acoustical analogy (see Tab. 1). To be compared to Fig. 15.

Figure 17 :

 17 Figure 17: A full vocal tract shown as a network of tracts. The small circles at the ends of the graph denote the control ports. The junctions are denoted by the circled variables η j . Here, tract T 2 and T 3 model the piriform fossae and tracts T 10 to T 12 model the nasal cavity. η 4 represents the nasal branching.

  into a DA-pH of form (3.3) using the incidence matrix DpGq of a directed graph G; and (ii) reformulates this fully implict differential-algebraic equation (DAE) into a semi-explicit system, of form Eq. (3.4), using a coordinate change[START_REF] Maschke | Port-Controlled Hamiltonian Systems: Modelling Origins and System Theoretic Properties[END_REF][START_REF] Van Der Schaft | Port-Hamiltonian Differential-Algebraic Systems[END_REF][START_REF] Wu | Port Hamiltonian System in Descriptor Form for Balanced Reduction: Application to a Nanotweezer[END_REF][START_REF] Cardoso | Port-hamiltonian Modeling and Control of a Fluid-structure System : Application to Sloshing Phenomena in a Moving Container Coupled to a Flexible Structure[END_REF].

  N s and j P I (3.6) where I " t1, . . . , N λ u " 1, N λ is the set of the vertex indexes. At last, we define a subset of control nodes N C Ă N (shown in pale red on Figs. 17 and 19), containing N C elements, and introduce N C -tuple C containing the indexes of all control nodes in N C such as C " ta P N : η a P N C u .

H

  k px k q ñ E " " e 1 , . . . , e k , . . . , e Ns ı . Concatenate the memoryless variables W " " w 1 . . . w Ns ı and associated laws ZpWq " " zpw 1 q . . . zpw Ns q ı .

  simplifying the interconnection -Accounting for the definition of P and the resulting split between constrained and unconstrained variables, we simplify Eq. (3.14) by observing that PC λ " " O pN λ ˆpnx´N λ qq I pN λ ˆNλ q ı

e θ pχ θ , χ λ q zpwq e λ pχ θ , χ λ q

Figure 22 :

 22 Figure 22: Directed graph of example 2: a tube of variable cross section composed of N t tracts.

Figure 23 :

 23 Figure 23: Spring-Damper system

FWFigure 24 :

 24 Figure 24: Equivalent electrical circuit of the shifted tract model augmented with simple resistors (see the extremities of the series branch). To be compared to Fig. 15.

Figure 25 :

 25 Figure 25: Two-tract system (N t " 2):

Figure 28 :

 28 Figure 28: Frequency responses of the two-tract model for both static phases: in solid blue, phase (1); in dashed black, phase (3).

Figure 29 :

 29 Figure 29: Frequency response of a uniform closed-open tube of total length 17 cm for N t " 1, 22 with an additional lumped radiation condition. The theoretical resonant frequencies (ideal closed-open resonator with an ideal radiation condition, see Ref. [CK13, 3.5.2, p. 273]) are denoted by the dashed vertical lines.

Figure 31 :

 31 Figure 31: Frequency response of the vocal tract for phases (1) and (3).

Figure 32 :

 32 Figure 32: Power flows of each tract during step (1) and (2).

Figure 33 :

 33 Figure 33: Flow rate at the lips and its time derivative for a simple coarticulation. Looking at the time derivative makes it possible to account for a simple radiation effect. Low frequency transients are present at t 1 and t 2 . The spectrogram of the bottom plot is shown on Fig. 34.

Figure 34 :

 34 Figure 34: Numerical results for the coarticulation [a]-to-[i]. From top to bottom: geometry at phase (1) (left plot) and at phase (3) (right plot); articulation coefficient, where 0 corresponds to phase (1) and 1 corresponds to phase (3); sum of the energy of the transverse flow of all tracts intensity normalised spectrogram; and energy per time frame. t 1 and t 2 are indicated by the dashed vertical lines.

Figure 35 :

 35 Figure 35: Input enthalpy and output mass rate at both extermities of the vocal tract.

Figure 36 :

 36 Figure 36: A generic three port junction between 3 systems. The cavity of the junction is shown in blue and its three connection boundaries are denoted with dashed red lines.

Figure 37 :

 37 Figure 37: The five different junction models derived in the chapter. Models (M1a-c) are identical w.r.t. the hypothesis: only the structure varies. Models (M2) and (M3) are obtained by applying linearising approximations.

  ¨n dS (5.7) where BΩ " S L Y S b Y S R Y S T with n being respectively ´ex ,´e y , e x and e y .

Figure 39 :

 39 Figure 39: Equivalent electrical electrical circuit of model (M1b). The current is analogous to the mass flow, and the voltage to the specific enthalpy. Inductors represent lumped fluid masses and the capacitor represents the compressible behaviour of the fluid.

Definition 2 (

 2 Linear time-invariant state and output equations): The linear time-invariant state-space representation of a dynamical system of state vector x reads $ & % dx dt " Axptq `Buptq yptq " Cxptq `Duptq (5.21) where A, B, C and D are constant matrices.

  R D S A P H M O D E L L I N G T O O L F O R F S I P R O B L E M S M O D E L L I N G I N F I N I T E D I M E N S I O N A L F S I P R O B L E M S W I T H C O L O R D I S T R I B U T I O N S I N L I G H T O F T H E P HF R A M E W O R K

Figure 40 :

 40 Figure 40: A simple fluid-structure interaction problem.

Definition 5 (

 5 Distribution,[START_REF] Schwartz | Méthodes mathématiques pour les sciences physiques[END_REF], Def. 2 p. 73): A continuous linear functional on the vector space D is called a distribution T .A distribution T P D 1 (dual space of D) assigns a function ϕ P D to a real number denoted by xT , ϕy where x¨, ¨y (or x¨, ¨yΩ ) is defined as the scalar product on D endowed with the following properties.Property 1 ([Sch65] (II, I; 8 and 9)):

Remark 15 (

 15 Empty intersection of two supports): If the support of ϕ P D and the support of T P D 1 have no points in common, then xT , ϕy " 0. differentiation of distributions [Sch65, p. 80] -We now define the spatial derivative of a distribution T .

¡Ω

  divpT ϕq dV ´xT , divpϕqy, where the first is null by the means of the divergence theorem[START_REF] Lipschutz | Vector Analysis and an Introduction to Tensor Analysis[END_REF] ¡ Ω divpT ϕq " ij BΩ T ϕ ¨n " 0, as ϕ P DpΩq is zero on BΩ.

∂ΩFigure 41 :

 41 Figure 41: The control domain Ω paved with subdomains ω i , i P 1, N ω

Figure 42 :

 42 Figure 42: The color function c i and its gradient along the associated subdomain ω i

ω π ρ ´gradp¨q 0

 0 v B and e B f " r π , ψs appear explicitely as inputs. D f corresponds to the mapping introduced in Ref. [DM13, Eq. (3.10)]. A dual operator D ‹ f exists, and is associated to output ports.

Figure 43 :

 43 Figure 43: Numerical results for the mass-spring-damper. The Python code is listed above.

Figure 44 :

 44 Figure 44: Numerical results for the adimensioned mass-spring-damper system. π and xs are dimensionless state variables, therefore no unit is indicated.

  that ρpx, y, tq " ρptq, we obtain E comp by multiplying U by the volume of the cavity E comp " P 0 V

  mr0s and hr0s are the initial values at t " 0 s. The associated outputs are ρ 0 ˙2 p h `hr0sq. (B.4.15) where ρp m, hq is a function that regenerates the fluctuations of the volumetric mass around its state of rest value ρ 0 . C A P P E N D I X O F C H A P T E R 4 c.1 the two-tract system: a toy model

Figure 45 :Figure 46 :Figure 47 :

 454647 Figure 45: Velocity input signal for the second tract

Figure 48 :¨´p 1

 481 Figure 48: Vowel synthesis: velocity signals applied to the walls.

  synthesis: velocity signals applied to the walls. . . . . . . . . . . . Correspondence table . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Table 2 Summary table for all junction models. . . . . . . . . . . . . . . . . . . . Table 3 Physical parameters used for the simulations of the two tract model . . . Table 4 Numerical parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

xvii L I S T O F TA B L E S

Table

1

  p0, tqe µ p0, tq ´ev pL, tqe µ pL, tq " f B e B

	leading to the effort variables
	δH δX	" ˜ρ0 Spxqv γP 0 µ {Spxqρ 2 0	¸" ˜ev e µ	¸, acoustic momentum linear density acoustic enthalpy
	and canonical pH formulation
						˜Bv	¸" «	ff	˜ev	¸(2.9)
							Bt				0	´gradp	¨q
							Bµ Bt			´divp ¨q	0	e µ
												J
	where ´divp ¨q and ´gradp ¨q are formally skew-adjoint according to Ref. [ZK14].
	The first line describes the conservation (specific [Dui+09, Sec. 3.4]) momentum,
	and the second line describes the conservation of mass.
	Finally, the power balance reads
	dHpX q dt	"	ż	Ω	δHpX q δX	BX Bt	dx	(2.9) "	ż	Ω	e J e dx
			(2.2) "	ż	Ω	´ev e µ	´divpe v q ¯˜´gradpe µ q	¸dx "	ż

Ω

´divpe v e µ q dx " e v

Table 1 :

 1 Correspondence table plicity and common usage in the voice modelling community. This also allows for easy comparisons with popular model of the lumped-parameter community [Mae82; BJK07; EL16]. It relies on a chosen analogy (see table

	Domains	Electronics	Solid mechanics	Fluid mechanics	Acoustics
	Quantities	voltage	force	specific total enthalpy	pressure
	(effort/flux)	current	velocity	mass rate	flow rate
	Component	inductor capacitor	mass spring	fluid mass compressibility	acoustic mass compliance
	Conservation laws	current law voltage law Newton's 2 nd law rigid connection	mass balance momentum balance

  Figure 26: Heights of T 1 and T 2 and total mass of air in the vocal tract versus time during the dynamical phase; t 1 and t 2 are shown in red dotted lines

		2.0	Height of tracts	Total mass	×10 -5
	Height (cm)	1.0 1.5		h 1 h 2		m 1 + m 2	2.5 3.0	Mass (kg)
			1	t 2	t 1	t 2
			Time (s)	Time (s)

  " rv mx , v c , v y , ρs x 2 " rν L , ν R , ν T , ms x 3 " r νL , νR , νT , ms x 4 " r νL , νR , νT , ms x 5 " rν L , ν R , ν T , ms

		(M1a)	(M1b)	(M1c)
	Section	5.3.1	5.3.2	5.3.3	5.3.4	5.3.5
	Hyp.
	App.				(App. 1)	(App. 1 & 2)
	Domain			Nonlinear acoustics		Linear acoustics
	Shifted
	Passive
	x x 1 Energy
	H k pxq

Table 2 :

 2 Summary table for all junction models.

  similar form is shown Ref. [OF04, Eq. (1.16)]. For higher order derivatives of T , see theorem 3 in Ref. [Sch65, p. 81]. examples of distribution -[Sch65, p.-77] The identity distribution 1 assigns a function ϕ to its integral on Ω, that is, x1 , ϕy " ş Ω ϕ. The Dirac distribution δ evaluates a function of D at r " 0 P Ω, which is denoted by xδ , ϕy " ϕp0q.

  The time derivative of c i pr, tq reads

			ż	T	x	Bc i Bt	, ϕy dt " rxc i , ϕys T 0	´żT	xc i ,	Bϕ Bt	y dt,
	where the first term vanishes, as ϕ is zero on BT . The second term, under Def. 6,
	reduces to		´żT	xc i ,	Bϕ Bt	y dt "	´żT	ˆx1 ,	Bϕ Bt	y ω i ˙dt,
	which, using the Leibniz-Reynolds transport theorem [Yam08, p. 26, Eq. (1.5.10)],
	gives (detailing the space-time variables)
	´żT	xc i p ' , tq ,	Bϕp ' , tq Bt	y ω i ptq dt "	´żT ¨d dt	¡
										ω i ptq

ϕpr, tq dV prq ‹ ' dt 0 `żT ¨ij Bω i ptq ϕpr, tq pv B pr, tq ¨npr, tqq dSpr, tq ‹ ' dt,

  B ρp ' , tq Bt `divp πp ' , tqq " ÿ iPF rDIV c i s πp ' , tq ´ÿ iPF rΓ c i πsp ' , tq (6.27)where we define the following distributions.Under Def. 8 and for c i pr, tq, we define the distribution Γ c i p that computes the flux of p P P induced by the motion of the boundary Bω i ptq, that isxrΓ c i psp ' , tq , ϕy " ij Bω i ptq c i p ' , tq pp ' , tqv B p ' , tqϕ ¨p´np ' , tqq dS, (6.28)where the integral over Bω i ptq, function of c i p, must be understood as the integral of the interior (inside ω i ) limit of p towards Bω i .For c i prq, we define the distribution DIV c i p P D 1 pΩ ˆT q that computes the flux of p through Bω

						(6.30)
	To do so we consider the scalar product x ρ , ϕy and examine its time derivative
	d dt	x ρp ' , tq , ϕy	"	dt d	x	ÿ
		0

i , xDIV c i p , ϕy " ij Bω i pc i pq ¨p´nq dS (6.29) where the integral over Bω i ptq assumes the same role as for Γ c i . Proof. The above theorem means that, for a test function ϕ P DpΩq, we have x B ρp ' , tq Bt `divp πp ' , tqq ´ÿ iPF rDIV c i s πp ' , tq `ÿ iPF rΓ c i πsp ' , tq , ϕy " 0. iPF ρi p ' , tqc i p ' , tq , ϕy " ÿ iPF d dt x ρi p ' , tqc i p ' , tq , ϕy,

  t. t is given by Eq. (6.17) on ω i , so that term 1

	becomes	1 "		ÿ iPF	¡ ω i ptq	B ρi Bt	ϕ dV "	iPF ÿ	¡ ω i ptq	´divp πi qϕ dV .
	We use the vectorial identity a b " a divpbq `gradpaq ¨b[GR07, p. 1051] to
	obtain
					¡
		1 "	ÿ				´divp πi ϕq `π i ¨gradpϕq dV
				iPF	ω i ptq
					¡				¡
		"	ÿ iPF	ω i ptq	´divp πi ϕq dV	`ÿ iPF	ω i ptq	πi ¨gradpϕq dV	.
									3	4
	The term 4 is given by
		4 "	ÿ
										¡
										Def. 9 "	π gradpϕq dV
										Ω
				¡
		"	Ω		divp πϕq dV	´¡ Ω	divp πqϕ dV
							5		6
	Term 5 is null as, under the divergence theorem and because ϕ is zero on BΩ
										ij
									5 "	πϕ ¨n dS " 0.	(6.31)
										BΩ
	Finally, given that d
	we have
		0 " 2 `3	`6
	d dt	x ρp ' , tq , ϕy "	ij Bω i ptq	ρ i v B ¨n ϕ dS	`ÿ iPF	¡ ω i ptq	´divp πi ϕq dV	´¡ Ω	divp πqϕ dV
	which, under the divergence theorem, gives back the quantity of the theorem.

iPF ¡ Ω c i πi gradpϕq dV as gradpϕq P DpΩq, " ¡ Ω ÿ iPF pc i πi q gradpϕq dV dt x ρ , ϕy " x B Bt ρ , ϕy as Ω is motionless and because ϕ P DpΩq,

  b.2 equations of dynamics b.2.1 Alternate Form of the Conservation of EquationThe classical form of the momentum conservation equation is is the entalpy and where we used the state equation p " ρ 2 B ρ pρq with pρq is the specific internal energy. Using these definitions and the identity v divpvq " div We consider a smooth function φpx, y, tq. Integrating over the domain Ω T ptq and using the Leibniz-Reynolds identity (also called the Reynold transport theorem) and the vectorial identity divpabq " gradpbq a `divpaqb yields This form is used to establish the dynamics of the transverse flow. We consider a test function with vectorial values. momentum conservation equation: alternate form This particular form is used to compute the dynamics of the transverse flow, as it directly uses the enthalpy. Its local formulation gives We now consider a smooth function with vectorial values. Using the same procedure than previously, the integral formulation reads:The state equation for pressure gives p " ρ 2 B ρ pρq. Averaging p over the domain of a tract yields In our case, we assume that the system (3.19) is under atmospheric conditions, meaning that the surrounding air is at rest at atmospheric pressure P 0 . Thus, the inputs reads y ‹ " Taking the equations of dynamics (3.19) where we substitute each energy variable by its state at rest value yields the following system of equations: Then, from Eq. (B.4.13c) we get m ‹ " S W h ‹ ρ 0 which also indicates ρ ‹ " m ‹ S W h ‹ " ρ 0 . Then, from Eq. (B.4.13d), one gets ν ‹ L " ´ν‹ R . Finally, using Eq. (B.4.13a) or Eq. (B.4.13b) yields ν ‹ L " ´ν‹ R " 0. The state at rest energy variables are:

					Bv Bt	`gradpψq " 0.
						Bρv B`d ivpρv b vq `gradppq "	(B.2.6)
	d dt	¡ Ω T ptq	For spatial domains with open boundaries, we use an alternative that includes the total specific entalpy ψ defined by ψ " 1 2 v ¨v `ωpρq with ωpρq " B Bρ pρ pρqq " pρq φv dV " ij BΩ φvpv B ¨nq dS `¡ Ω T ptq Bφ Bt ¨v dV ´¡ Ω T ptq divpφψq dV divpφqψ dV . `¡ Ω T ptq `ρ B pρq Bρ " pρq `p{ρ (B.2.7) (B.2.11)
	b.3 weak formulation of the pressure
								´v2	ūnder
								2
				the irrotational hypothesis, the conservation of momentum equation is written
			such as This form can be seen in [vM02]. Bv Ω T ptq Bt ă p ą" 1 V ¡ p dV " 1 `gradpψq " 0. Ω T ptq ¡ V ρ 2 B ρ pρq dV .	(B.3.12)	(B.2.8)
	b.4 shifted version of the tract model
				"	ı		"	ı
	ψ ‹ L ψ ‹ R F ‹ W compression energy (see App. B.1).	"	Bρ Bt 0 0 ´P0 S W " ´divpρvq. under the choices for the
							"	ı
	We assume that a state at rest x ‹ T "	ν ‹ L ν ‹ R Π ‹ y m ‹ h ‹	exists such as
			d dt	¡ Ω T ptq d dt ν ‹ L " 0 " ´em ´vy φρ dV " ij BΩ pφρq pv B ´vq ¨n dS Π ‹ y m ‹ d dt ν ‹ R " 0 " e m `vy Π ‹ y m ‹	`¡ Ω T ptq	ρ	Bφ Bt	dV	gradpφq ρv dV (B.4.13a) `¡ Ω T ptq (B.2.9) (B.4.13b)
			d dt momentum conservation equation The integral formulation reads: Π ‹ y " 0 " d dt pm ‹ q Π ‹ y m ‹ `F ‹ W ´eh " F ‹ W ´eh d dt m ‹ " 0 " e ν L `eν R	(B.4.13c) (B.4.13d)
	d dt First, Eq. (B.4.13e) gives that Π ‹ ¡ Ω T ptq φρv dV " ij BΩ φρvpv B ¨nq dS d h ‹ " 0 " e Πy . dt y " 0.	`¡ Ω T ptq	Bφ Bt	¨ρv dV	´¡ Ω T ptq	divpρv b vqφ dV (B.4.13e)	(B.2.10)
								´ij BΩ	pφ ¨n dS	`¡ Ω T ptq	divpφqp dV .

b.2.2 Integral Formulation with a Time-Varying Domain mass conservation equation The mass conservation equation reads: d dt x ‹ T " 0.

Table 3 :

 3 Physical parameters used for the simulations of the two tract model

	Parameter	Value
	Sampling frequency F s 1 ¨10 5 Hz
	Machine's precision	2 ¨10 ´16
	DG	1.5 ¨10 ´8
	Table 4: Numerical parameters

Please refer to section

2.2 in Ref.[vJ14] for further details on Dirac structures in pHs.

See Ref.[vJ14, p. 167] for more details.

I would like to thank Denis Matignon for the fruitful exchanges on this example.

For numerical purposes, the condition δxn ‰ 0 can be replaced by |δxn| ě where is the numerical precision of the chosen floating point representation. A possible refinement is proposed in this chapter.

Available at https://github.com/pyphs/pyphs.

Where the numerical precision is e.g. » 2 ¨10 ´16 for double floats.

Subscripts L and R stand for the left and right open surfaces, respectively; the W subscript stands for wall, the top surface.

Also called static pressure or thermodynamic pressure.

i.e., an interconnection matrix only filled with values 0 and ˘1.

These form functions correspond to the linear functions associated to the left and right surfaces, respectively, which can be thought as the nodes of a mesh.

laminar resistance of a laminar flow within a circular duct[START_REF] Kenneth | Airflow and Turbulence Noise for Fricative and Stop Consonants: Static Considerations[END_REF].

Please refer to Ref.[START_REF] Wetzel | Power-balanced Modelling of the Vocal Tract: A Recast of the Classical Lumped-parameter Model[END_REF] for an example that involves a non-canonical interconnexion matrix. Be aware that this reference provides a particular case and that further computations could be needed for another pHs.

They will remain unchanged after the coupling procedure, like the external port pv W , F W q of the wall in the case of the tract model shown on Eq. (3.19)

A first experimental study on the nonlinear behaviour of this model has been conducted in Ref.[START_REF] Dalmont | Experimental determination of the equivalent circuit of an open side hole: linear and non linear behaviour[END_REF]

Quantities are expressed as per unit of mass.

In this appendix, we give the volumetric internal energy density U pρq " ρ pρq.

This is only valid in the absence of cavitation effects, as there are no holes between the fluid and solid phases

See the immersed boundary methods[START_REF] Griffith | Immersed Methods for Fluid-Structure Interaction[END_REF] or the "artificial implicit compressibility added coupling" method[START_REF] Olivier | Numerical Strategies for Partitioned Fluid-structure Interaction Algorithms[END_REF]. See Ref[START_REF] Landajuela | Coupling Schemes for the Fsi Forward Prediction Challenge: Comparative Study and Validation[END_REF] for a comparison of coupling methods.

Inspired by Refs. [GTR05; BSW15].

& & &

beautiful village of Frauenchiemsee (Germany). /?lang=en https://hamecmopsys.ens2m.fr. https://anr.fr/Projet-ANR-15-CE23-0007. https://websites.isae-supaero.fr/infidhem/.

Similarily as in the first experiment, we consider three temporal phases:

(1) a static phase t P rt 0 " 0 s, t 1 " 0.07 ss, where the vocal tract articulates the vowel [a];

(2) a dynamical phase t P rt 1 , t 2 " 0.08 ss, where the walls of the vocal tract are driven to reach the geometry of vowel [i];

(3) a static phase t P rt 2 " 0.08 s, t max " 0.7 ss where the vocal tract articulates the vowel [i].

The data for the vocal tract geometries of each vowel have been taken from Ref. [START_REF] Arnela | MRI-Based Vocal Tract Representations for the Three-Dimensional Finite Element Synthesis of Diphthongs[END_REF].

For each tract, the driving velocity signal at the wall is given by the boxcar function, the amplitude of which is ratio of the motion distance to the time allotted for the movement (0.01 s in this simulation). To avoid harsh discontinuities at the start and at the end of the movement, we smoothed out the extremeties of the velocity step with sinusoidal ramps (see Fig. 47 and Fig. 48 in App. C).

For both static phases, we plot the simulated geometry on Fig. 30. [i]

Figure 30: Simulated geometry of the vocal tract for N t " 15. In blue, the geometry of vowel [a] associated with phase (1). In orange, the geometry of the vowel [i] associated with phase (3). During the dynamics phase (2), the geometry changes linearily with time.

To compute the associated frequency responses (see Fig. 31), we feed an entahlpy impulse at the glottis at t imp,1 " 0.035 s and at t imp,2 " 0.62 s.

Step 3 : the kinematic (ii) and dynamic (iii) conditions are specified, leading to the fully coupled formulation of the FSI problem.

An important feature is that the physical hypotheses and coupling conditions are treated separately, making this methodology modular. It is easy to differentiate the coupling hypotheses, acting at the boundaries, and the physical hypotheses, governing the dynamics of each phase within their subdomain. On top of this feature, the resulting models are physically interpretable, and can be discretised with a pHcompatible spatial discretisation method [START_REF] Luiz | A Partitioned Finite Element Method for power-preserving discretization of open systems of conservation laws[END_REF] to produce power-balanced finitedimensional models.

However, this methodology is still limited in its applications, as it only accounts for one-dimensional problems composed of two phases.

To address these issues, we attempt, in this last chapter, to generalise the work of Diagne and Maschke [START_REF] Diagne | Port Hamiltonian Formulation of a System of Two Conservation Laws with a Moving Interface[END_REF] in the three-dimensional case and for an arbitrary number of phases. We choose a distribution theory setting, where we consider weak forms of the physical variables and equations with precise definitions.

In order to make this preliminary work fit in a chapter, we simplify some elements. In Ref. [START_REF] Diagne | Port Hamiltonian Formulation of a System of Two Conservation Laws with a Moving Interface[END_REF], the considered dynamical systems are abstract, as they may be any one-dimensional system of two conservation equations. This choice enables the authors to focus on mathematical matters, namely the investigation of the formal pH structure of the coupled dynamical system. In this preliminary work, we consider a physical problem with a simple infinite-dimensional description of a perfect fluid, from which we develop our methodology. This allows for a concise and physical take on the subject, in view of inspiring forthcoming studies to propose new and more intuitive formulations.

The chapter is organised as follows. In section 6.1, we describe the principle of levelset methods, then cast their main tool, the indicator function (or color function), as a distribution. In section 6.2, we show how to obtain an uncoupled formulation of the dynamics of the fluid phases, augmented with color distributions. In section 6.3, we give a few perspectives about a possible port-Hamiltonian formulation of this preliminary work.

color distributions

Indicator Functions in Level-Set Methods

Level-set methods [MOS92; Set99] are used in image processing, computer vision, and CFD for FSI problems or multiphase flow (see Ref. [START_REF] Osher | The Level Set Methods and Dynamic Implicit Surfaces[END_REF] for more details). This class of methods relies on the use of an indicator function, also called level-set function or color function, to enforce the geometric conditions [START_REF] Richter | Fluid-structure Interactions[END_REF] and describe the moving boundaries enclosing each phase. We use the most basic definition of color function, detailed hereafter. In the following code, we implement the mass-spring-dashpot example as a Core object by specifying the energy variables (lines 47 and 51), the associated Hamiltonians (lines 48 and 52), algebraic components (lines 55-56) and external ports (line 59). The full model is then assembled and the interconnexion matrix is specified.

The rest of the code is related to the simulation and the plot of the numerical results. ---------------------------------------------------------------# --------------------------------------------------------------- # Adding the components to the pHs model core . add_storages ( pi , hm ) core . add_storages (x , hk ) core . add_dissipations ( v_m , z ) core . add_ports (F , v ) ---------------------------------------------------------------# --------------------------------------------------------------- ---------------------------------------------------------------# --------------------------------------------------------------- # ----------------------------------------------------------------# --------------------------------------------------------------- ---------------------------------------------------------------# --------------------------------------------------------------- where Spxq is divided as

where we ommit the dependence of the block matrices on x for clarity.

Consider a set of adimensionalising coefficients rx s1 , . . . , x snx s and the associated dimensionalising matrix X s " diag px s1 , . . . , x snx q of size pn x ˆnx q (A.2.3) so that the state vector can be written as