Passive modelling of the electrodynamic loudspeaker: from the Thiele–Small model to nonlinear port-Hamiltonian systems

The electrodynamic loudspeaker couples mechanical, magnetic, electric and thermodynamic phenomena. The Thiele/Small (TS) model provides a low frequency approximation, combining passive linear (multiphysical or electric-equivalent) components. This is commonly used by manufacturers as a reference to specify basic parameters and characteristic transfer functions. This paper presents more refined nonlinear models of electric, magnetic and mechanical phenomena, for which fundamental properties such as passivity and causality are guaranteed. More precisely, multiphysical models of the driver are formulated in the core class of port-Hamiltonian systems (PHS), which satisfies a power balance decomposed into conservative, dissipative and source parts. First, the TS model is reformulated as a linear PHS. Then, refinements are introduced, stepby-step, benefiting from the component-based approach allowed by the PHS formalism. Guaranteed-passive simulations are proposed, based on a numerical scheme that preserves the power balance. Numerical experiments that qualitatively comply with measured behaviors available in the literature are presented throughout the paper.


Introduction
The electrodynamic loudspeaker is a non ideal transducer. Its dynamics is governed by intricate multiphysical phenomena (mechanical, magnetic, electric and thermodynamic), a part of which involves nonlinearities responsible for audio distortions [1][2][3]. As a first example, the viscoelastic properties of the suspension material result in long-term memory (linear) and hardening spring effect (nonlinear). Second, the ferromagnetic properties of the solid iron core in the voice-coil result in eddy current losses (linear) and magnetic saturation (nonlinear). Third, the resistance of the coil wire converts a part of electrical power into heat. This modifies material properties and, eventually, can cause irreversible damages. Such phenomena must be modeled and considered in the design of real-time distortion compensation [4][5][6][7] and that of burn-out protection [2,8].
In the sequel, we refer to this set of parameters as the Thiele-Small model. It combines passive linear models of elementary physical components (see Fig. 1) and provides a low-frequency linear time-invariant approximation for low-amplitude excitation on short period. This (multiphysical or electric-equivalent) parametric model is commonly used by manufacturers as a reference to specify basic parameters and characteristic transfer functions.
Various refinements of this reference model have been proposed, both in the frequency domain and the time domain [3,[13][14][15][16][17][18]. In particular, the lumped-parameter approach [1,2,19,20] consists in modeling the dependence of Thiele-Small parameters on some selected physical quantities (e.g. position-dependent stiffness). However, fundamental physical properties are usually not guaranteed by the proposed mathematical models, and their physical interpretation is not always obvious. Examples are model causality in the context of frequency domain simulation based on Fourier transform and model passivity in the context of arbitrary polynomial approximation of constitutive relations of materials. Obviously, this is also the case for gray-box modelling based on Volterra and Wiener/Volterra series [21][22][23][24] or nonlinear ARMAX [25]).
To circumvent those difficulties, we introduce the port-Hamiltonian systems formalism [26][27][28] as a systematic framework for the nonlinear modeling, simulation and control of loudspeakers. Port-Hamiltonian systems are state-space representations that satisfy a power balance structured into conservative, dissipative and external (/source) parts. This structure can be described by acausal graphs (such as electronic circuits or classical electromechanical analogies) which allow modular exploration or local model refinements, while guarantying passivity. Such modular processes are facilitated by existing methods [29] that automatically generate the governing equations associated with a given multiphysical system described by a network of predefined or user-defined components. It is then possible to run guaranteed-passive simulations based on a numerical method that preserves the power balance and its structure in the discrete-time domain, from which passivity and stability properties stem. The objective of this work is to model well-known multiphysical phenomena occurring in loudspeakers as a set of port-Hamiltonian structures and components in view of system identification and correction.
The elementary phenomena that are considered in this paper are concerned with mechanical, magnetic and electric phenomena, as well as their coupling. They are known to be responsible for significant audio distortions (see e.g. [1,2] and references therein): The materials used for the suspension (S) exhibit, first, a combination of the behaviors of elastic solids and viscous fluids ( [30], Sect. 1.2), inducing long time shape memory (creep effect, see e.g. Fig. 1 from [30] and Fig. 11 from [31]) and second, nonlinear stressstrain characteristics so that the restoring force is not proportional to the elongation [1,18], with maximal instantaneous excursion q sat that corresponds to the breakdown of the material. The electromechanical coupling (back e.m.f. and Lorentz force on the diaphragm D) depends on the coil (C) position and on the magnetic flux in the pole piece (P). The former is modulated by the movement of the coil which acts as an electromagnet that modifies the latter. The materials used for the pole piece P exhibit nonlinear magnetic excitation-induction curve so that the equivalent current in the coil is not proportional to its magnetic flux. Also, a maximal magnetic flux / sat is reached (flux saturation), corresponding to the global alignment of the microscopic magnetic moments (see [32] Sect. 1 and Ref. [33]).
Most of the magnetic materials (iron, cobalt, etc.) possess high electric conductivity. The application of a variable magnetic induction induces currents, namely eddy-currents, in a plane orthogonal to the field lines (see [34], Sect. 1.1.2). This has three effects: a power is dissipated due to the natural resistivity of the material (Joule effect), eddy-currents induces their own magnetic field (added inductive effect), and they oppose to the original induction (Lens's law), which pushes the field lines toward the boundary (magnetic skin effect). This paper is structured as follows. Section 2 recalls the Thiele-Small model (model 0) which recast as a port-Hamiltonian system after a short introduction to the formalism. This model 0 serves as a basis to elaborate two refined loudspeaker models (model 1 and model 2). Section 3 focuses on refinements of the mechanical part (model 1). In particular, a passive-guaranteed nonlinear model based on the Kelvin-Voigt description of viscoelastic material associated with the suspension (S) is provided. Section 4 focuses on refinements of the electromagnetic part (model 2). In particular, we provide a lumped parameter non linear description of the full magnetic circuit (M, P) coupled with the electronic circuit and the mechanical system. These passive models can straightforwardly be combined to describe all these refinements in a single model. Simulations are presented throughout the paper, focusing on the effect of each phenomena, separately. In practice, they are all 1 produced by the PyPHS software [35].

The Thiele-Small model revisited in the port-Hamiltonian formalism
This section is devoted to the construction of the base model (model 0) that is progressively refined in the remaining of the paper. First, an overview of the functioning of the electrodynamic loudspeaker is presented and the standard Thiele-Small modeling is recalled. Second, the port-Hamiltonian framework is recalled. Third, the Thiele-Small model is recast as a port-Hamiltonian system (model 0). Finally, time-domain simulations are performed and results are compared in the frequency-domain to transfer functions expected from filter theory.

Standard Thiele/Small model
The basic functioning of a boxed loudspeaker such as the one depicted in Figure 1 is as follows. A voice-coil (C) submitted to an input voltage (I) is immersed in a magnetic field imposed by a permanent magnet (M) in the air gap (G) of a magnetic path (pole piece P), so that the coil (C) is subjected to the Lorentz force. The coil (C) is glued to a large diaphragm (D) which is maintained by a flexible suspension (S). The diaphragm (D) is responsible for the coupling with the acoustical field (A). 1 Simulation code are available here: https://afalaize.github.io/ posts/loudspeaker1/ The standard description of the dynamics of this system is referred as the Thiele-Small modeling, introduced in the early seventies [9][10][11][12]. The electrical part (C) includes the electrical resistance of the coil wire R C and the linear approximation of the coil behavior with inductance L C . The mechanical part (C, D, S, A) is modeled as a damped harmonic oscillator with mass M CDA (coil, diaphragm and additional mass due to acoustic radiation), linear approximation of the spring effect K SA (suspension and additional stiffness due to air compression in the enclosure) and fluidlike damping with coefficient R SA (frictions and acoustic power radiation). The magnetic part (M, P, G, C) reduces to a constant force factor B'.
The corresponding set of ordinary differential equations are derived by applying Kirchhoff's laws to the electrical part (C) and Newton's second law to the mechanical part (D, S, A): with v I the input voltage, i C the coil current and q D the diaphragm's displacement from equilibrium. The electromechanical coupling terms are the back electromotive force

Port-Hamiltonian formalism
The port-Hamiltonian (pH) formalism introduced in the 90's [26] is a modular framework for the passive-guaranteed modeling of open dynamical systems. In this paper, we consider the following class formulated as a multiphysical, component-based, differential algebraic state-space representation (as in [29]). Definition 2.1 (Port-Hamiltonian Systems, PHS). The class of PHS under consideration is that of differential algebraic state-space representations with input u 2 R ny , state x 2 R nx , output y 2 R ny , that are structured according to energy flows and described by (see [29] for details and [26][27][28] for more general formulations of PHS ): where w 2 R nw stands for dissipation variables with dissipation law zðwÞ 2 R nw , HðxÞ 2 R þ is the energy storage function (or Hamiltonian) with gradient ðrHðxÞÞ i ¼ oH ox i , K 2 R nxÂnw , G x 2 R nxÂny , G w 2 R nwÂny , and where (iii) J x 2 R nxÂnx , J w 2 R nwÂnw and J y 2 R nyÂny are skewsymmetric matrices, so that J T = ÀJ.
System (3) proves passive for the outgoing power P S = u T y according to the following power balance: This proves the passivity and hence the asymptotic stability of (3) in the sense of Lyapunov ( [36], Sect. 4).
Remark 2.2 (Energy storage). The Hamiltonian considered in this work does not depend explicitly on time so that it coincides with the total energy in the system:

The Thiele-Small model as a PHS
The Thiele-Small modeling from Section 2.1 can be regarded as the interconnection of a resistance-inductance circuit with a mass-spring-damper system, through a gyrator that describes the reversible energy transfer from the electrical domain to the mechanical domain as depicted in Figure 2 and detailed in Section B.2.

Description
This system includes n x = 3 storage components (inductance L C , mass M CDA and stiffness K SA ), n w = 2 dissipative components (electrical resistance R C and mechanical damping R SA ) and n y = 1 port (electrical input v I ). The state

Port-Hamiltonian formulation
The above PHS quantities are related with quantities in the Thiele-Small model (1) and (2) as follows: Electrical part: dx 1 dt This yields the following port-Hamiltonian reformulation of (1)-(2): The associated port-Hamiltonian structure (3) is given in Table 1.

Simulation results
Simulations are performed following the passive-guaranteed numerical method associated with the pH structure (3) and recalled in Appendix A. Time domain simulations are shown in Figure 3. Transfer function computed from time with s the complex frequency is shown in Figure 4. Notice the (numerical) power balance is satisfied. The model 0 in Table 1 is refined in the sequel to cope with the phenomena listed in the introduction.

Refined mechanics
In this section, the model 0 from Section 2.3 is refined to cope with creep effect and nonlinear stress-strain relation attached to the suspension material (S). First, we detail the modeling of the creep effect based on Kelvin-Voigt model of viscoelastic material [30,31]. This results in a linear PHS. Second, the hardening suspension effect is included. This results in a nonlinear PHS (model 1). Third, simulation results are shown.

Suspension creep
The creep effect is a long-term memory effect due to the shape memory of the suspension material (see e.g. Fig. 1 from [30] and Fig. 11 from [31]) and heat relaxation of the fluid in the enclosure (see e.g. [17]).
In this work, we consider the standard Kelvin-Voigt model for viscoelastic materials. The resulting (linear) mechanical subsystem is depicted in Figure 5 and is recast in this subsection as a port-Hamiltonian system (3). Note the procedure given below allows to formulate other creep models (e.g. from the literature given above) as port-Hamiltonian systems.

Description of the creep model
Viscoelastic materials exhibit combination of elastic solids behaviors and viscous fluid behaviors. Let R be the coefficient of viscosity for a damper (N Á s Á m À1 ) and K the modulus of elasticity for a spring element (N Á m À1 ) with characteristic frequency x ¼ K R (Hz) and associated creep Rs where s is the complex Laplace variable ðReðsÞ > 0Þ, and where q(s), f K (s) and f R (s) are the Laplace transforms of the elongation and the two restoring forces, respectively.
The Kelvin-Voigt modeling of the creep effect is constructed by connecting a linear spring with same stiffness K in parallel with a damper R (see [37] and Sect. 4 from [38]). The elongation is the same for both elements q kv = q and forces sum up f kv = f K + f R . The corresponding compliance is The modeling of materials that exhibits several relaxation times s n ¼ 2p x n is achieved by chaining N Kelvin-Voig elements (see [20,38,39] and Fig. 1 from [37]). Each element contributes to the total elongation q chain ¼ P N n¼1 q n , and every elements experience the same force Here, we consider three elements to restore (i) a primary instantaneous response to a step force with stiffness K 0 and (ii) a long time viscoelastic memory with characteristic time s ve ¼ 2p x 1 . The compliance of this viscoelastic model is The parameters are K 0 , K 1 and s ve with . A possible strategy to tune jointly the fK i g i¼1;2 is to introduce a dimensionless parameter P K 2 ð0; 1Þ to partition the Thiele-Small stiffness K SA between K 0 and K 1 : Table 1. Port-Hamiltonian formulation (3) for the Thiele-Small structure as depicted in Figure 2, with magnetic flux in the coil / C , diaphragm position q D and momentum p M ¼ M CDA dq D dt . The physical parameters are given in Tables D1.   Storage  State: Energy: Dissipation Variable: Law: Note this choice ensure that T ve ðsÞj sve¼0 ¼ K SA =s for all 0 < P K < 1 (i.e. the combination of K 0 and K 1 restores the Thiele-Small stiffness K SA if the creep effect is neglected s ve = 0).

Port-Hamiltonian formulation
The creep model (7) corresponds to the parallel connection of (i) a linear spring K 0 and (ii) a linear spring K 1 serially connected to a dashpot R 1 (see Fig. 5). This mechanical subsystem includes n x = 3 storage components (mass M CDA , primary stiffness K 0 , secondary stiffness K 1 ), n w = 2 dissipative components (damper R SA , secondary damper R 1 ), and n y = 1 port (Lorentz force f L ). The state dt , and the primary and secondary elongations q 0 and q 1 (respectively). The Hamiltonian is the sum of (i) the kinetic energy H M ðx 1 Þ ¼ 2M CDA , and (ii) the primary and secondary potential energies For these definitions, the interconnection in Figure 5 yields:  Table 1.
Physical parameters are given in Table D1. The input voltage v I is a 100 Hz sine wave with increasing amplitude between 0 and 50 V. The sampling rate is F s = 96 kHz.  Table 1 (phs) and computed from filter theory (target). Physical parameters are given in Table D1. The input voltage v I is a 1 V white noise and the sampling rate is F s = 96 kHz. Figure 5. Small-signal modeling of the mechanical part which includes: the total mass M CDA (diaphragm, coil and additional mass due to acoustic radiation), the fluid-like damper R SA (mechanical friction and small signal approximation for the acoustic power radiation), primary stiffness K 0 and Kelvin-Voigt modeling of the creep effect (K 1 , R 1 ), with diaphragm position q D , primary elongation q 0 and creep elongation q 1 . Parameters are given in Tables D1 and D2.
This system is recast as a port-Hamiltonian system (3) for the structure in Table 2 and the parameters in Table D2.

Suspension hardening and model 1
For large displacement, the suspension behaves like an hardening spring (see e.g. [18,40]). This should occur for instantaneous displacements, so that only the primary stiffness K 0 in Table 2 is affected. First, the mechanical subsystem from previous section is changed to cope with this phenomenon. Second, the resulting nonlinear mechanical part is included in loudspeaker model 0 to build the loudspeaker model 1.

Model description
The primary stiffness K 0 in Table 2 is modified into a nonlinear spring that exhibits a phenomenological saturation for an instantaneous elongation q 0 = ±q sat (symmetric). The associated constitutive law (C1)-(C3) in Appendix C is It yields the restoring force f 0 (q 0 ) = K 0 c SA (q 0 ) for the initial stiffness K 0 . It corresponds to the addition of a saturating term that does not contribute around the origin, thus preserving the meaning of parameter K 0 (small signal behavior). The associated storage function (C4) and (C5) is Note that any storage function suitable to a particular material could be used (such as for example those given in Appendix C), and that the modular structure of the port-Hamiltonian formalism allows to change the storage function without modifying the interconnection matrix.

Port-Hamiltonian system and Model 1
The port-Hamiltonian formulation of the loudspeaker model 1 includes creep effect and hardening suspension. It is obtained by (i) replacing the potential energy K 0 Table 2 by the nonlinear storage function (11) and (ii) connecting the mechanical port f L to the RL circuit describing the electromagnetic part as in Section 2.3. This results in the structure given in Table 3 with parameters  in Tables D1 and D2.  (3) for the proposed small signal model of the mechanical part in Figure 5 driven by the Lorentz force f L , with diaphragm position q D , momentum p M ¼ M CDA dq D dt , primary elongation q 0 , and creep elongation q 1 . Parameters are given in Tables D1 and D2, with Storage State: Energy: Output:  (3) for the model 1 depicted in Figure 6. The linear stiffness K SA is replaced by the Kelvin-Voigt modeling of the creep effect from Section 3.1 in serial connection with the nonlinear spring described in Section 3.2, with diaphragm position q D , momentum p M ¼ M CDA dq D dt , primary elongation q 0 , and creep elongation q 1 . The nonlinear potential energy H SA sat ðq 0 Þ is given in (11). Parameters are given in Tables 4 and D1, with Storage State: Energy: Law: Output:

Simulation results
The numerical method used to simulate the mechanical subsystem (Tab. 2) and the model 1 (Tab. 3) is detailed in Appendix A. The results obtained for physical parameters in Tables D1 and D2 are commented below.

Creep effect
It is expected that the viscoelastic behavior of the suspension material results in a frequency-dependent compliance, i.e. the suspension at low frequencies must appear softer than for the Thiele/Small prediction (see e.g. [17], Fig. 12). Model 1 allows the recovery of this effect as shown in Figure 7. The corresponding long time memory depicted in Figure 8 is in accordance with measurements in e.g. Figure 1 from [30] and Figure 11 from [31].

Nonlinear suspension
The hardening effect associated with the nonlinear stress-strain characteristic of the suspension material is clearly visible in Figure 9 where the primary elongation is reduced for higher value of the shape parameter P S sat . This reduces the total displacement q D and momentum p M ¼ M CDA dq D dt , while the creep elongation is almost unchanged.

Refined electromagnetic
In this section, the model 0 from Section 2.3 is refined to account for effects of flux modulation, electromagnetic coupling, ferromagnetic saturation and eddy-current losses attached to the electromagnetic part (voice-coil C, magnet M, ferromagnetic path P and air gap G). First, the proposed modeling is described. Second, this model is recast as a port-Hamiltonian system. Third, simulation results are presented.

Model description
The classical lumped elements modeling of loudspeakers electrical impedance includes the electrical DC resistance of the wire R C serially connected to a non-standard inductive effect, referred as lossy-inductor. The simplest refinement of the Thiele/Small modeling is the so-called LR-2 model, which uses a series inductor connected to a second inductor shunted by a resistor. Several refinements to this model are available in the litterature [16,[41][42][43][44].
The modeling of the loudspeakers electrical impedance proposed in this work is depicted in Figure 10. The coil winding acts as an electromagnetic transducer (gyrator) that realizes a coupling between the electrical and the magnetic domains, according to the gyrator-capacitor approach (see Appendix B.3 and Refs. [45,46]). The electrical domain includes the linear resistance R C of the coil wire (same as Thiele-Small model) and a constant linear inductance associated with the leakage magnetic flux that does not penetrate the pole piece (P). The flux in the magnetic   Table 2: Compliance in the frequency domain (diaphragm displacement in response to the Lorentz   path is common to (i) a nonlinear magnetic capacitor associated with energy storage in air gap (G, linear) and ferromagnetic (P, nonlinear), (ii) a linear magnetic dissipation associated with eddy-currents losses in the path (P) and (iii) a constant source of magnetomotive force associated with the permanent magnet (Ampere model).

Coil model Leakage inductance
A single leakage flux / leak = S leak b leak independent of the position q D is assumed for every of the N C wire turns (see Fig. 10a), with S leak the annular surface between the coil winding and the ferromagnetic core, computed as where 0 < a leak < 1 is the fraction of the coil section not occupied by the magnetic core. According to (B1), the linear magnetic capacity of the air path is C leak ¼ S leak l 0 ð1þn air Þ 2 A C with A C the height of the coil wire turns, l 0 the magnetic permeability of vacuum and n air the magnetic susceptibility of air. From (B6), this corresponds to an electrical inductance with state x leak = N C / leak and storage function H leak ðx leak Þ ¼ x leak 2 2L leak , for the inductance L leak ¼ N 2 C C leak . We define the characteristic frequency

Electromagnetic coupling modulation
The electromagnetic coupling between the coil (C) and the path (P) depends on the he number n P of wire turns effectively surrounding the pole piece. For small negative excursions q D < 0 every wire turns participate to the coupling (n P ' N C ) and for large positive excursions the coil leaves the pole piece (n P ' 0). We propose a phenomenological sigmoid relation n P : q D 7 !n P ðq D Þ: with n P (q À ) ' 90% N C and n P (q + ) ' 10% N C (see Fig. 11). Figure 9. Simulation of the model 1 in Table 3 depicted in Figure 6, for the parameters in Tables D1 and D2 (except P S sat indicated in the legend). The input voltage is a 100 Hz sine wave with increasing amplitude between 0V and 50V. The sampling rate f s = 96 kHz. The power balance is shown for P S sat ¼ 10 only. Notice q D = q 0 + q 1 . Figure 10. Proposed modeling of the electromagnetic circuit, which includes: the coil wire resistance R C , the linear inductance associated with the leakage flux / leak , the electromagnetic transduction with n P the number of wire turns around the magnetic path, the magnetic energy storage in the ferromagnetic path described by the nonlinear induction-excitation curve w PG (/ PG ) from (15), the linear dissipation associated with eddycurrents in the pole piece, and the constant source of magnetomotive force w M due to the magnet from (17).

Ferromagnetic saturation Nonlinear storage
The storage of magnetic energy in the magnetic circuit is spread over the pole piece (P) and the air gap (G). Assuming no leakage flux in the pole piece, those elements are crossed by the same magnetic flux / PG (see Fig. 10a and Refs. [45,46]). The corresponding averaged inductions are with S P the average section crossed by the magnetic flux in the pole piece and S G the section of the flux in the air gap (see Fig. 10a). This corresponds to the serial connection of two magnetic capacitors: the first one is associated with the air gap G with linear constitutive law (as for the leakage flux / leak ); the second one is associated with the pole piece P and cannot be described by a linear magnetic capacity due to the magnetic saturation that occurs in ferromagnetic material (see [32], Sect. 1). Those two serially-connected magnetic capacitors can merge into a single nonlinear capacitor that restores the total magnetomotive force w PG (/ PG ). In this work, we consider the tangent-like constitutive relation detailed in Appendix C with flux saturation / sat = S P b sat , where b sat depends on the specific magnetic material. From (C1) to (C3), the constitutive law w PG (/ PG ) = c PG (/ PG ) is given by where the coefficient P PG lin includes the contributions of both air and pole piece material, and P PG sat is a function shape parameter that depends on the specific material used for the pole piece. The associated (positive definite) storage function (C4) and (C5) is given by Again, any storage function suitable to a particular magnetic material could be used (such as those given in Appendix C), and the modular structure of the port-Hamiltonian formalism allows to change the storage function without modifying the interconnection matrix.

Steady state behavior
The permanent magnet is modeled as a constant source of magnetomotive force w M (Ampere model [46]). This drives the magnetic flux in the path to an equilibrium (steady-state, ss) / PG = / ss for which the magnetomotive force exactly compensates that of the magnet: The associated steady-state magnetic capacity is the inverse of the linear approximation of w PG (/ PG ) at / ss : with C ss ¼ L P n P ð0Þ 2 and L P = L C ÀL leak .

Eddy-currents losses
Besides the magnetic saturation, the pole piece is affected by the combination of capacitive and resistive effects due to eddy-currents, resulting in frequencydependent losses. This phenomenon is well described by a linear fractional order magnetic capacity (see [3,34,[47][48][49][50] and [48], part 5). This is out the scope of the present work Figure 11. Plot of the position-dependent effective number of wire turns n P (q D ) involved in the electromagnetic coupling from (13) with q + = 5 mm and q À in the legend.  Table 4. Elements common to model 0 in Figure 2 are shaded. The coil inductance L C is replaced by the electromagnetic circuit from Figure 10b which includes the leakage inductance L leak and the magnetic path (pole piece P and air gap G). and is postponed to a follow-up paper. Here, we consider a magnetic resistance R ec (X À1 ) with magnetic impedance Since we consider a single magnetic flux in the pole piece / ec = / PG , this impedance is serially connected to the magnetic capacity described in Section 4.1.2. The resulting structure is depicted in Figure 10b. Defining x P = (R ec C ss ) À1 (Hz) and s P ¼ 2p x P (s), the resulting electrical impedance T P ðsÞ ¼ t P i C is

Blocked electrical impedance
The current i C is common to (i) the resistor R C , (ii) the leakage inductance L leak , and (iii) the impedance associated with the magnetic path in the coil core T P ðsÞ. For a blocked , the total steady-state electrical impedance T C ðsÞ ¼ v I ðsÞ i C ðsÞ measured at the coil terminals is given by The DC value (s = 0) is given by the resistance R C . In the high frequency range, the impedance is governed by the leakage inductance T C ði xÞ $ bracket is the contribution of the proposed magnetic circuit.

Position-dependent force factor
The gyrator that restores the Lorentz force f L with corresponding back electromotive force v L is given by (see with coil velocity v C ¼ dq D dt and ' C the length of coil wire effectively subjected to the magnetic field B. The apparent magnetic field B is modulated by the coil displacement (neglected in this work), and the length ' C depends on the coil position (see Figs. 2.5-2.8 in [2] and Fig. 5 in [1]). We propose a parametric plateau function ' C : q D 7 !' C ðq D Þ to model the latter phenomenon: where ' C 0 is the total length of the coil, Q ' describes the overhang of the coil with respect to the magnetic path (see Fig. 13a; and Sect. 3.1.2 from [1]), and P ' is a shape parameter (see Fig. 13b).

Port-Hamiltonian formulation
The proposed loudspeaker modeling that includes electromagnetic phenomena (model 2) corresponds to the replacement of the inductance L C in model 0 by the electromagnetic circuit described in previous section (compare Figs. 2 and 12). It includes (i) the resistance-inductance circuit R C À L leak serially connected to (ii) the magnetic circuit associated with the core of the coil and the magnet. This involves n x = 4 storage components (inductance L leak , capacity c PG , mass M CDA , and stiffness K SA ), n w = 3 dissipative components (resistances R C , R SA and R ec ), and n y = 2 ports (voltage v I and magnetomotive force w M ). The state is x ¼ ðx leak ; / PG ; p M ; q D Þ T with the state associated with leakage flux x leak ¼ N C / leak , and the

Hamiltonian is HðxÞ
and H PG sat given in (16). The dissipation variable is w ¼ i C ; dq D dt ; w PG À Á T with linear dissipation law zðwÞ ¼ diagðR C ; R SA ; R ec À1 Þ w. According to (14), the magnetic induction in the air gap involved in the electromechanical coupling (B3) is b G ¼ / PG S G . The length of wire effectively subjected to the induction field is ' C ðx 4 Þ given in (24). The number of wire turns effectively surrounding the pole piece involved in the electromagnetic coupling is n P (x 4 ) given in (13). With these definitions, the port-Hamiltonian formulation (3) of the loudspeaker model with the refined electromagnetic part (model 2 depicted in Fig. 10) is given in Table 4.

Simulation results
The numerical method used to simulate model 2 is detailed in Appendix A. Physical parameters are given in Tables D1 and D3. In each case, the initial condition is the steady-state / PG (t = 0) = x 2 (t = 0) = / ss . Figure 13. Effective length of coil wire ' C subjected to the magnetic field B as defined in (24), with coil position q D and total wire length ' C 0 ¼ 10 m. Notice P ' ¼ 0 corresponds to ' C ¼ ' C 0 , 8q D 2 R which restores the linear case. (a) Effect of the overhang parameter Q ' with P ' ¼ 10. (b) Effect of the shape parameter P ' with Q ' ¼ 5 mm.

Eddy-current losses
The effect of the characteristic time s ec = 2pR ec C ss due to eddy-currents in the pole piece is illustrated by imposing several DC input voltages v I , here À50 V, and +50 V (see evolution of flux / PG in Fig. 14), with the coil blocked at q D = 0 m and C ss kept fixed, so that only R ec varies with s ec . In each case, the magnetic flux in the path is driven to a new steady state / PG 6 ¼ / ss . The long-term effects and the influence of the characteristic time s PG are clearly visible.

Core saturation
The evolution of the small signal impedance with the steady-state is shown in Figure 15. First, the DC input voltages are imposed during 0.5 s. Second, a small signal noise is applied to evaluate the new steady-state blockedimpedance. We see an evolution in the high-frequency response according to the transfer function in (22). The associated value for is given in C ss from (18) and the other parameters are given in Table D3.

Position-dependent electromagnetic coupling
To illustrate the effect of coil position on the electrical impedance, position q D in model 2 (Tab. 4) is fixed to À1 cm (inside), 0 cm (equilibrium) and +1 cm (outside). Due to the position-dependent effective number of coil wire (13), this changes the inductance according to (B6). Results are shown in Figure 16, in accordance with measurements in e.g. Fig. 6 from [1].

Flux-dependent force factor
The force factor in model 2 B' ¼ / PG S G ' C ðq D Þ is modulated by the coil position (same as model 0) and the magnetic flux Table 4. Blocks associated with the port-Hamiltonian formulation (3) for the loudspeaker model 2 depicted in Figure 10, where the Lorentz force factor is B'ðxÞ ¼ x2 S G ' C ðx 4 Þ with the magnetic induction in the air gap / PG =S G and the positiondependent effective wire length ' C ðq D Þ defined in (24). See definitions and notations in Section 4. Storage State: Energy: x ¼ Figure 14. Simulation of the loudspeaker model 2 in Table 4: Normalized flux / PG / ss in response to a ±50 V step voltage. The initial flux is / PG (t = 0) = / ss and the coil is blocked at q D = 0 m. The sample-rate is 96 kHz. Figure 15. Simulation of the loudspeaker model 2 in Table 4: in the pole piece P and air gap G. This is clearly visible in the results of Figure 17: we observe that the force factor can be larger than predicted by the Thiele/Small modeling. Notice that the power balance is fulfilled.

Conclusion
In this paper, a set of structures and components have been proposed to model well-known multiphysical phenomena occurring in loudspeakers, in view of system identification and correction. In particular, a finite-dimensional, power-balanced and passive-guaranteed time-domain formulation of viscoelastic and eddy-currents phenomena (linear) and material properties (stress-strain and b-h characteristics, nonlinear) have been derived. Those models are given in the framework of port-Hamiltonian systems, which decomposes the system into conservative, dissipative and source parts. The numerical method used for the simulations preserves this decomposition and thus is unconditionally stable. Numerical results that qualitatively comply with measured behaviors available in the literature have been presented.
The two loudspeaker models 1 and 2 have been developed independently of each other. This permits to illustrate the particular effect of each phenomenon on the loudspeaker dynamics. Now, their interconnection to form a global, multiphysical model that copes with all the phenomena covered in this work is straightforward, due to the modular nature of the port-Hamiltonian systems.
The first perspective of this work is to achieve DSP simulation-based real-time audio distortion compensation, based on the preliminary work in [7]. This requires the development of a parameter estimation method dedicated to the port-Hamiltonian structure. A second perspective is to include other phenomena that have not been considered here, such that the fractional dynamics associated with viscoelastic materials and eddy-currents, the acoustical load and the thermal evolution of the system. For all these issues, the modular structure of the proposed port-Hamiltonian models could be further exploited.   Table 4: Variation of impedance when the coil is blocked in different positions, hence changing the number of coil wire turns around the path n P (q D ) and the inductance according to (37). The flux is initially at / PG (t = 0) = / ss . The sample-rate is 96 kHz. Figure 17. Simulation of the model 2 in Table 4 depicted in Figure 12, for the parameters in Tables D1 and D3. The input voltage is a 100 Hz sine wave with increasing amplitude between 0 V and 50 V. The sampling rate is 96 kHz. The power balance is shown for the model 2 only. The force factor B corresponds to the product of the induction in air gap b G from (14) with position-dependent effective length from (24). discrete time-domain to preserve passivity. This is achieved by numerical schemes that provide a discrete version of the chain rule for computing the derivative of E ¼ H x. This is the case of Euler scheme, for which first order approximation of the differential applications dxðt; dtÞ ¼ dx dt ðtÞ dt and dH(x, dx) = rH(x) T dx on the sample grid t ¼ k dt; k 2 Z are given by dxðk; dtÞ ¼ xðk þ 1Þ À xðkÞ;

Acknowledgments
ðA1Þ For mono-variate storage components (H ðxÞ ¼ P n x n¼1 H n ðx n Þ), the solution can be built element-wise with the n-th coordinate given by For pH systems composed of a collection of linear energy storing components with quadratic Hamiltonian H n ðx n Þ ¼ which restores the midpoint rule. For nonlinear case, (A3) does not coincide with the mid-point rule anymore, still preserving passivity due to equation (A6) which does not assume the system linearity.
Remark A.1 (Nonlinear solver). The proposed method guarantees the passivity of the discrete time model provided the resulting nonlinear implicit equations are solved exactly. This depends on the nonlinear solver at hand but not on the proposed method.
In this section, we give the elements for the modeling of lumped electromagnetic systems in the pH formalism. First, the closed-form expression for energy storage is recalled. Second and third, we give the port-Hamiltonian formulation of electromechanical and electromagnetic coupling, respectively.

Definitions
The magnetic phenomena are described by two complementary fields, namely, the applied magnetic excitation h and the induced magnetic flux density b(h), which is somewhat the response of a given material to a given excitation. The induction b is defined as the superposition of the magnetization of vacuum j 0 (h) and the magnetization of matter j(h) due to microscopic magnetic moments attached to the atoms of the body (see Eq. (1.6) from [32] and Eq. (6) from [33] where we neglect the magnetization of vacuum so that h(b) = j À1 (b). The magnetic induction flux / is defined as the flux of the magnetic induction field through a given surface S : is the magnitude of the field b(t) that we assume constant over S and normal to the cross section. The magnetomotive force w is defined as the circulation of h along a closed b-field line C with length ' C : wðbðtÞÞ ¼ H C hðbðtÞÞ d' ¼ ' C Á hðbðtÞÞ, where we assume h(b(t)) constant along C.

Energy storage
The variation of magnetic energy density (locally) stored in a sample of magnetic material is dE dt ¼ hðbÞ db dt (see e.g. [33] for details). Again assuming that h is constant over the b-field line and b constant over a cross section of the material, the variation of the total energy for a sample with length ' C and cross section S is dE dt ¼ S ' C hðbÞ db dt . Rewriting the preceding relation for the magnetic flux / = Sb and the magnetomotive force wðbÞ ¼ ' C hðbÞ yields dE dt ¼ w / S À Á d/ dt . Thus, we can select / as the state associated with the storage of magnetic energy with E(t) = H mag (/(t)) and we identify H 0 mag ð/Þ ¼ w / S À Á : with the total energy variation d dt H mag ¼ w d/ dt .

B.2 Electromechanical coupling
Consider several windings of a conductive wire with section S W , total length ' W , position q W and velocity vector v W = v W e W with constant direction e W and magnitude v W ¼ dq W dt . This conductor is immersed in a magnetic induction field b with constant direction orthogonal to e W and constant magnitude B. The current is i W ¼ RR S W q q v q dS for the electric charge density q q moving inside the wire at velocity v q = v q e q with unitary vector e q normal to a cross section of the wire. A wire element with length d' is subjected to the Lorentz force df L ¼ This force is orthogonal to the velocity v q + v W so that the associated mechanical power is dP L ¼ df L Á ðv q þ v W Þ ¼ 0. Integrating along the wire, one gets defining the Lorentz force f L and the back electromotive force (voltage) v L . Notice the transfer is reversible and conservative in the sense that the outflow of energy from the electrical domain P elec ¼ i W v L equals the inflow of the mechanical domain P meca ¼ v W f L . This corresponds to a gyrator with ratio B ' W : since the interconnection matrix is skew-symmetric.

B.3 Electromagnetic coupling: the gyrator-capacitor approach
The gyrator-capacitor approach introduced in the late sixties [45,46] is an easy way to develop electronic analog of magnetic circuits. It has been considered in [14] for the modeling of the loudspeaker. In this approach, a coil is divided in a gyrator (wire turns) and a magnetic energy storage (coil core).
The dynamics of a magnetic field can be described by two complementary macroscopic quantities: the magnetic induction flux / and the magnetomotive force (mmf) w (see Sect. B.1). The electromagnetic transfer for a single wire turn stands from (i) Faraday's law of electromagnetic induction that relates the electromotive force (voltage t) to the variation of the magnetic flux in the wire turn t ¼ d/ dt ; and (ii) Ampere's theorem that relates the mmf to the current in the wire w = i (see [45,46]). Considering a coil (C) with N C wire turns around the path (P), these relations restore a gyrator with ratio N C : Denoting by s 2 C the Laplace variable, the correspondence between an impedance seen in the electrical domain Z elec ðsÞ ¼ t C ðsÞ i C ðsÞ and its counterpart in the magnetic domain Z mag ðsÞ ¼ w C ðsÞ s / C ðsÞ ¼ N 2 C i C ðsÞ t C ðsÞ is given by As a result, if the path (P) is modeled by a magnetic capacity C P , e.g. from the linearization of H 0 mag ð/Þ in (B1), the equivalent electrical inductance is L C ¼ N 2 C C P . Notice the interconnection (B5) is conservative: P elec = P mag with P elec = v C i C the power outgoing the electrical domain and P mag ¼ d/ C dt w C the power incoming the magnetic domain.