© A. Falaize and T. Hélie, Published by EDP Sciences, 2020

1 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–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–7] and that of burn-out protection [2, 8].

The basic reference set of parameters describing the electrodynamic loudspeaker is that of Thiele–Small [9–12]. 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–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–24] or nonlinear ARMAX [25]).

To circumvent those difficulties, we introduce the port-Hamiltonian systems formalism [26–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 stress–strain 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¹ produced by the PyPHS software [35].

2 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.

2.1 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).

The standard description of the dynamics of this system is referred as the Thiele–Small modeling, introduced in the early seventies [9–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 fluid-like damping with coefficient R_SA (frictions and acoustic power radiation). The magnetic part (M, P, G, C) reduces to a constant force factor $B\mathcal{l}$.

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):

$$v_{\mathrm{I}}\left(t\right)=v_{\mathcal{L}}\left(t\right)+R_{\mathrm{C}} i_{\mathrm{C}}\left(t\right)+L_{\mathrm{C}}\frac{\mathrm{d}{i}_{\mathrm{C}}\left(t\right)}{\mathrm{d}t},$$(1)

$$M_{\mathrm{CDA}}\frac{{\mathrm{d}}^2{q}_{\mathrm{D}}\left(t\right)}{\mathrm{d}{t}^2}=f_{\mathcal{L}}\left(t\right)-R_{\mathrm{SA}}\frac{\mathrm{d}{q}_{\mathrm{D}}\left(t\right)}{\mathrm{d}t}-K_{\mathrm{SA}} q_{\mathrm{D}}\left(t\right),$$(2)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 (voltage) $v_{\mathcal{L}}=B\mathcal{l} \frac{\mathrm{d}{q}_{\mathrm{D}}}{\mathrm{d}t}$ and the Lorentz force $f_{\mathcal{L}}=B\mathcal{l} i_{\mathrm{C}}$.

2.2 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 $\mathbf{u}\in {\mathbb{R}}^{n_{\mathit{y}}}$ , state $\mathit{x}\in {\mathbb{R}}^{n_{\mathit{x}}}$ , output $\mathbf{y}\in {\mathbb{R}}^{n_{\mathit{y}}}$ , that are structured according to energy flows and described by (see [29] for details and [26–28] for more general formulations of PHS):

$$\left(\begin{array}{l}\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}\\ \mathbf{w}\\ \mathbf{y}\end{array}\right)=\underset{\mathbf{J}}{\underbrace{\left(\begin{array}{ccc}{\mathbf{J}}_{\mathbf{x}}& -\mathbf{K}& {\mathbf{G}}_{\mathbf{x}}\\ {\mathbf{K}}^{\mathrm{T}}& {\mathbf{J}}_{\mathbf{w}}& {\mathbf{G}}_{\mathbf{w}}\\ -{\mathbf{G}}_{\mathbf{x}}^{\mathrm{T}}& -{\mathbf{G}}_{\mathbf{w}}^{\mathrm{T}}& {\mathbf{J}}_{\mathbf{y}}\end{array}\right)}}\left(\begin{array}{l}\nabla \mathrm{H}\left(\mathbf{x}\right)\\ \mathbf{z}\left(\mathbf{w}\right)\\ \mathbf{u}\end{array}\right),$$(3) where $\mathbf{w}\in {\mathbb{R}}^{n_{\mathbf{w}}}$ stands for dissipation variables with dissipation law $\mathbf{z}\left(\mathbf{w}\right)\in {\mathbb{R}}^{n_{\mathbf{w}}}$ , $\mathrm{H}\left(\mathbf{x}\right)\in {\mathbb{R}}_{+}$ is the energy storage function (or Hamiltonian) with gradient $\left(\nabla \mathrm{H}\right(\mathbf{x}){)}_i=\frac{\mathrm{\partial }\mathrm{H}}{\mathrm{\partial }{x}_i}$ , $\mathbf{K}\in {\mathbb{R}}^{n_{\mathbf{x}}\times n_{\mathbf{w}}}$ , ${\mathbf{G}}_{\mathbf{x}}\in {\mathbb{R}}^{n_{\mathbf{x}}\times n_{\mathbf{y}}}$ , ${\mathbf{G}}_{\mathbf{w}}\in {\mathbb{R}}^{n_{\mathbf{w}}\times n_{\mathbf{y}}}$ , and where

1. the storage function H(x) is positive semidefinite H(x) ≥ 0 with H(0) = 0 and positive definite Hessian matrix $\left[{\mathcal{H}}_{\mathrm{H}}\right(\mathbf{x}){]}_{i,j}=\frac{{\mathrm{\partial }}^2\mathrm{H}}{\mathrm{\partial }{x}_i,\mathrm{\partial }{x}_j}(\mathbf{x})$ (see examples inAppendix C),

2. the dissipation lawz(w) is null at originz(0) = 0 with positive definite Jacobian matrix $\left[{\mathcal{J}}_{\mathbf{z}}\right(\mathbf{w}){]}_{i,j}=\frac{\mathrm{\partial }{z}_i}{\mathrm{\partial }{w}_j}(\mathbf{w})$ , implying that the dissipated power is P_D (w) = z(w)^Tw ≥ 0, P_D (0) = 0,

3. ${\mathbf{J}}_{\mathbf{x}}\in {\mathbb{R}}^{n_{\mathbf{x}}\times n_{\mathbf{x}}}$ , ${\mathbf{J}}_{\mathbf{w}}\in {\mathbb{R}}^{n_{\mathbf{w}}\times n_{\mathbf{w}}}$ and ${\mathbf{J}}_{\mathbf{y}}\in {\mathbb{R}}^{n_{\mathbf{y}}\times n_{\mathbf{y}}}$ are skew-symmetric matrices, so thatJ^T = −J.

System (3) proves passive for the outgoing power P_S = u ^T y according to the following power balance:

$${\left(\begin{array}{l}\nabla \mathrm{H}\left(\mathbf{x}\right)\\ \mathbf{z}\left(\mathbf{w}\right)\\ \mathbf{u}\end{array}\right)}^{\mathrm{T}}\left(\begin{array}{l}\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}\\ \mathbf{w}\\ \mathbf{y}\end{array}\right)=\frac{\mathrm{dH}}{\mathrm{d}t}\left(\mathbf{x}\right)+\underset{\ge 0}{\underbrace{{\mathrm{P}}_{\mathrm{D}}\left(\mathbf{w}\right)}}+{\mathrm{P}}_{\mathrm{S}}=0.$$(4)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: $E=\mathrm{H}\circ \mathbf{x}:t\mapsto H\left(\mathbf{x}\right(t\left)\right)=E\left(t\right)\in {\mathbb{R}}_{+}.$

2.3 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 $\mathbf{x}={({\varphi }_{\mathrm{C}},p_{\mathrm{M}},q_{\mathrm{D}})}^{\mathrm{T}}$ consists of the magnetic flux in the coil ϕ_C, mass momentum $p_{\mathrm{M}}=M_{\mathrm{CDA}} \frac{\mathrm{d}{q}_{\mathrm{D}}}{\mathrm{d}t}$ and diaphragm position q_D. The Hamiltonian is the sum of the electrodynamic energy ${\mathrm{H}}_{\mathrm{L}}\left(x_1\right)=\frac{x_1^2}{2L_{\mathrm{C}}}$, the kinetic energy ${\mathrm{H}}_{\mathrm{M}}\left(x_2\right)=\frac{x_2^2}{2M_{\mathrm{CDA}}}$ and the potential energy ${\mathrm{H}}_{\mathrm{K}}\left(x_3\right)=K_{\mathrm{SA}}\frac{x_3^2}{2}$. The dissipation variable is $\mathbf{w}={\left(i_{\mathrm{C}},\frac{\mathrm{d}{q}_{\mathrm{D}}}{\mathrm{d}t}\right)}^{\mathrm{T}}$ with linear dissipation law $\mathbf{z}\left(\mathbf{w}\right)=\mathrm{diag}(R_{\mathrm{C}},R_{\mathrm{SA}})\mathbf{w}$.

Port-Hamiltonian formulation

The above PHS quantities are related with quantities in the Thiele–Small model (1) and (2) as follows:

• Electrical part: $\frac{\mathrm{d}{x}_1}{\mathrm{d}t}=L_C\frac{\mathrm{d}{i}_C\left(t\right)}{\mathrm{d}t}, \hspace{1em}{H}_L^{\text{'}}\left(x_1\right)=i_C\hspace{1em}\mathrm{and}\hspace{1em}{z}_1\left(w_1\right)=R_C{i}_C;$

• Mechanical part: $\frac{\mathrm{d}{x}_2}{\mathrm{d}t}=M_{\mathrm{CDA}}\frac{{\mathrm{d}}^2{q}_D\left(t\right)}{\mathrm{d}{t}^2}, \hspace{1em}{H}_M^{\text{'}}\left(x_2\right)=\frac{\mathrm{d}{q}_D}{\mathrm{d}t},\hspace{1em} H_K^{\text{'}}\left(x_3\right)=K_{\mathrm{SA}} q_D\hspace{1em}\mathrm{and}\hspace{1em}{z}_2\left(w_2\right)=R_{\mathrm{SA}}\frac{\mathrm{d}{q}_D}{\mathrm{d}t};$

with the back electromotive force $B\mathcal{l} {\mathrm{H}}_{\mathrm{M}}^{\mathrm{\prime }}\left(x_2\right)=v_{\mathcal{L}}$ and the Lorentz force $B\mathcal{l} {\mathrm{H}}_{\mathrm{L}}^{\mathrm{\prime }}\left(x_1\right)=f_{\mathcal{L}}$. This yields the following port-Hamiltonian reformulation of (1)–(2):

$$\begin{array}{c}\frac{\mathrm{d}{x}_1}{\mathrm{d}t}=-B\mathcal{l} \frac{\mathrm{\partial }\mathrm{H}}{\mathrm{\partial }{x}_2}\left(x_2\right)-z_1\left(w_1\right)+u_1,\\ \frac{\mathrm{d}{x}_2}{\mathrm{d}t}=B\mathcal{l} \frac{\mathrm{\partial }\mathrm{H}}{\mathrm{\partial }{x}_1}\left(x_1\right)-\frac{\mathrm{\partial }\mathrm{H}}{\mathrm{\partial }{x}_3}\left(x_3\right)-z_2\left(w_2\right),\end{array}$$(5)with coil velocity $v_{\mathrm{C}}=\frac{\mathrm{\partial }\mathrm{H}}{\mathrm{\partial }{x}_2}\left(x_2\right)=\frac{\mathrm{d}{x}_3}{\mathrm{d}t}=w_2$ and current $i_{\mathrm{C}}=\frac{\mathrm{\partial }\mathrm{H}}{\mathrm{\partial }{x}_2}\left(x_2\right)=w_1=-y_1$. 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 domain simulation $\mathcal{T}\left(s\right)=\mathrm{abs}\left(\frac{v_{\mathrm{C}}\left(s\right)}{v_{\mathrm{I}}\left(s\right)}\right)$ 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.

3 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.

3.1 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.

3.1.1 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⁻¹) and K the modulus of elasticity for a spring element (N · m⁻¹) with characteristic frequency $\omega =\frac{K}{R}$ (Hz) and associated creep time $\tau =\frac{2\pi }{\omega }$ (s). Their respective compliance in the Laplace domain are ${\mathcal{T}}_K=\frac{q\left(s\right)}{f_K\left(s\right)}=\frac{1}{K}$ and ${\mathcal{T}}_R=\frac{q\left(s\right)}{f_R\left(s\right)}=\frac{1}{\mathrm{Rs}}$ where $s$ is the complex Laplace variable $\left(\mathfrak{R}e\right(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

$$\begin{array}{ccc}{\mathcal{T}}_{\mathrm{kv}}\left(s\right)& =& \frac{q_{\mathrm{kv}}\left(s\right)}{f_{\mathrm{kv}}\left(s\right)}={\left(K \left(1+\frac{s}{\omega }\right)\right)}^{-1}.\end{array}$$(6)

The modeling of materials that exhibits several relaxation times ${\tau }_n=\frac{2\pi }{{\omega }_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_{\mathrm{chain}}={\sum }_{n=1}^N\mathrm{}{q}_n$, and every elements experience the same force $f_{\mathrm{chain}}=f_1=\cdots =f_N$. Here, we consider three elements to restore (i) a primary instantaneous response to a step force with stiffness K₀ and (ii) a long time viscoelastic memory with characteristic time ${\tau }_{\mathrm{ve}}=\frac{2\pi }{{\omega }_1}$. The compliance of this viscoelastic model is

$${\mathcal{T}}_{\mathrm{ve}}\left(s\right)=\frac{1}{K_0}+{\left(K_1 \left(1+\frac{s}{{\omega }_1}\right)\right)}^{-1}, {\omega }_1=\frac{K_1}{R_1}.$$(7)

The parameters are K₀, K₁ and τ_ve $\left(\mathrm{with} R_1=\frac{K_1{\tau }_{\mathrm{ve}}}{2\pi }\right)$. A possible strategy to tune jointly the $\{K_i{\}}_{i=\mathrm{1,2}}$ is to introduce a dimensionless parameter $P_K\in \left(\mathrm{0,1}\right)$ to partition the Thiele–Small stiffness K_SA between K₀ and K₁:

$$K_0=\frac{K_{\mathrm{SA}}}{1-P_K}, \hspace{1em}{K}_1=\frac{K_{\mathrm{SA}}}{P_K}.$$(8)Note this choice ensure that ${{\mathcal{T}}_{\mathrm{ve}}\left(s\right)|}_{{\tau }_{\mathrm{ve}}=0}=K_{\mathrm{SA}}/s$ for all 0 < P _K  < 1 (i.e. the combination of K₀ and K₁ restores the Thiele–Small stiffness K_SA if the creep effect is neglected τ_ve = 0).

3.1.2 Port-Hamiltonian formulation

The creep model (7) corresponds to the parallel connection of (i) a linear spring K₀ and (ii) a linear spring K₁ serially connected to a dashpot R₁ (see Fig. 5). This mechanical subsystem includes n _x  = 3 storage components (mass M_CDA, primary stiffness K₀, secondary stiffness K₁), n _w  = 2 dissipative components (damper R_SA, secondary damper R₁), and n _y  = 1 port (Lorentz force f_L). The state $\mathbf{x}=(p_{\mathrm{M}},q_0,q_1{)}^{\mathrm{T}}$ includes the mass momentum $p_M=M_{\mathrm{CDA}}\frac{\mathrm{d}{q}_{\mathrm{D}}}{\mathrm{d}t}$, and the primary and secondary elongations q₀ and q₁ (respectively). The Hamiltonian is the sum of (i) the kinetic energy ${\mathrm{H}}_{\mathrm{M}}\left(x_1\right)=\frac{x_1^2}{2M_{\mathrm{CDA}}}$, and (ii) the primary and secondary potential energies ${\mathrm{H}}_0\left(x_2\right)=K_0\frac{x_2^2}{2}$ and ${\mathrm{H}}_1\left(x_3\right)=K_1\frac{x_3^2}{2}$ (respectively). The dissipation variable is $\mathbf{w}={\left(\frac{\mathrm{d}{q}_{\mathrm{D}}}{\mathrm{d}t},f_{R_1}\right)}^{\mathrm{T}}$ with linear dissipation law $\mathbf{z}\left(\mathbf{w}\right)=\mathrm{diag}(R_{\mathrm{SA}},R_1^{-1})\cdot \mathbf{w}$. The input/output are $\mathbf{u}=(f_{\mathcal{L}}{)}^{\mathrm{T}}$ and $\mathbf{y}={\left(\frac{\mathrm{d}{q}_{\mathrm{D}}}{\mathrm{d}t}\right)}^{\mathrm{T}}$. For these definitions, the interconnection in Figure 5 yields:

$$\begin{array}{c}\frac{\mathrm{d}{x}_1}{\mathrm{d}t}=-\frac{\mathrm{\partial }\mathrm{H}}{\mathrm{\partial }{x}_2}\left(x_2\right)-z_1\left(w_1\right)+u_1,\\ \frac{\mathrm{d}{x}_2}{\mathrm{d}t}=\frac{\mathrm{\partial }\mathrm{H}}{\mathrm{\partial }{x}_1}\left(x_1\right)-z_2\left(w_2\right),\\ \frac{\mathrm{d}{x}_3}{\mathrm{d}t}=z_2\left(w_2\right).\end{array}$$(9)This system is recast as a port-Hamiltonian system (3) for the structure in Table 2 and the parameters in Table D2.

3.2 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₀ 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.

3.2.1 Model description

The primary stiffness K₀ in Table 2 is modified into a nonlinear spring that exhibits a phenomenological saturation for an instantaneous elongation q₀ = ±q_sat (symmetric). The associated constitutive law (C1)–(C3) in Appendix C is

$$c_{\mathrm{SA}}\left(q_0\right)=q_0+\frac{4 P_{\mathrm{sat}}^{\mathrm{S}}}{4-\pi } \left(\tan \left(\frac{\pi .q_0}{2 q_{\mathrm{sat}}}\right)-\frac{\pi .q_0}{2q_{\mathrm{sat}}}\right).$$(10)

It yields the restoring force f₀(q₀) = K₀ c_SA(q₀) for the initial stiffness K₀. It corresponds to the addition of a saturating term that does not contribute around the origin, thus preserving the meaning of parameter K₀ (small signal behavior). The associated storage function (C4) and (C5) is

$${\mathrm{H}}_{\mathrm{sat}}^{\mathrm{SA}}\left(q_0\right)=K_0\left(\frac{q_0^2}{2}-\frac{8P_{\mathrm{sat}}^{\mathrm{S}} q_{\mathrm{sat}}}{\pi (4-\pi )} \left(\ln \left|\cos \left(\frac{\pi q_0}{2 q_{\mathrm{sat}}}\right)\right|+\frac{1}{2}{\left(\frac{\pi q_0}{2q_{\mathrm{sat}}}\right)}^2\right)\right)\ge 0.$$(11)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.

3.2.2 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\frac{q_0^2}{2}$ in Table 2 by the nonlinear storage function (11) and (ii) connecting the mechanical port $f_{\mathcal{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.3 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_{\mathrm{sat}}^{\mathrm{S}}$. This reduces the total displacement q_D and momentum $p_{\mathrm{M}}=M_{\mathrm{CDA}}\frac{\mathrm{d}{q}_{\mathrm{D}}}{\mathrm{d}t}$, while the creep elongation is almost unchanged.

4 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.

4.1 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–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 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).

4.1.1 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

$$S_{\mathrm{leak}}=\frac{\pi D_{\mathrm{C}}^2}{4} (2 {\alpha }_{\mathrm{leak}}-{\alpha }_{\mathrm{leak}}^2),$$(12)where 0 < α_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_{\mathrm{leak}}=\frac{S_{\mathrm{leak}} {\mu }_0 (1+{\xi }_{\mathrm{air}})}{2 A_{\mathrm{C}}}$ with A_C the height of the coil wire turns, μ₀ the magnetic permeability of vacuum and ${\xi }_{\mathrm{air}}$ the magnetic susceptibility of air. From (B6), this corresponds to an electrical inductance with state x_leak = N_C ϕ_leak and storage function ${\mathrm{H}}_{\mathrm{leak}}\left(x_{\mathrm{leak}}\right)=\frac{{x_{\mathrm{leak}}}^2}{2L_{\mathrm{leak}}}$, for the inductance $L_{\mathrm{leak}}=N_{\mathrm{C}}^2 C_{\mathrm{leak}}$. We define the characteristic frequency ${\omega }_{\mathrm{C}}=\frac{R_{\mathrm{C}}}{L_{\mathrm{leak}}}$ (Hz).

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_{\mathrm{P}}:q_{\mathrm{D}}\mapsto n_{\mathrm{P}}\left(q_{\mathrm{D}}\right)$:

$$n_{\mathrm{P}}\left(q_{\mathrm{D}}\right)=N_{\mathrm{C}} {\left(1+\exp \left(\frac{4 q_{\mathrm{D}}-2 (q_{-}+q_{+})}{q_{+}-q_{-}}\right)\right)}^{-1},$$(13)with n_P(q₋) ≃ 90% N_C and n_P(q₊) ≃ 10% N_C (see Fig. 11).

4.1.2 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

$$\begin{array}{c}{b}_{\mathrm{P}}=\frac{{\varphi }_{\mathrm{PG}}}{S_{\mathrm{P}}},\\ b_{\mathrm{G}}=\frac{{\varphi }_{\mathrm{PG}}}{S_{\mathrm{G}}},\end{array}$$(14)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 ψ_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 ψ_PG(ϕ_PG) = c_PG(ϕ_PG) is given by

$$\begin{array}{l}{c}_{\mathrm{PG}}\left({\varphi }_{\mathrm{PG}}\right)=P_{\mathrm{lin}}^{\mathrm{PG}}\left({\varphi }_{\mathrm{PG}}+\frac{4 P_{\mathrm{sat}}^{\mathrm{PG}}}{4-\pi } \left(\tan \left(\frac{\pi {\varphi }_{\mathrm{PG}}}{2 {\varphi }_{\mathrm{sat}}}\right)-\frac{\pi {\varphi }_{\mathrm{PG}}}{2{\varphi }_{\mathrm{sat}}}\right)\right),\end{array}$$(15)where the coefficient $P_{\mathrm{lin}}^{\mathrm{PG}}$ includes the contributions of both air and pole piece material, and $P_{\mathrm{sat}}^{\mathrm{PG}}$ 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

$${\mathrm{H}}_{\mathrm{sat}}^{\mathrm{PG}}\left({\varphi }_{\mathrm{PG}}\right)=P_{\mathrm{lin}}^{\mathrm{PG}}\left(\frac{{\varphi }_{\mathrm{PG}}^2}{2}-\frac{8P_{\mathrm{sat}}^{\mathrm{PG}} {\varphi }_{\mathrm{sat}}}{\pi (4-\pi )}\left(\ln \left|\cos \left(\frac{\pi {\varphi }_{\mathrm{PG}}}{2 {\varphi }_{\mathrm{sat}}}\right)\right|+\frac{1}{2}{\left(\frac{\pi {\varphi }_{\mathrm{PG}}}{2{\varphi }_{\mathrm{sat}}}\right)}^2\right)\right).$$(16)

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 ψ_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:

$${\psi }_{\mathrm{PG}}\left({\varphi }_{\mathrm{ss}}\right)=-{\psi }_{\mathrm{M}}.$$(17)

The associated steady-state magnetic capacity is the inverse of the linear approximation of ψ_PG(ϕ_PG) at ϕ_ss:

$$C_{\mathrm{ss}}={\left({\frac{\mathrm{\partial }{c}_{\mathrm{PG}}}{\mathrm{\partial }{\varphi }_{\mathrm{PG}}}|}_{{\varphi }_{\mathrm{PG}}={\varphi }_{\mathrm{ss}}}\right)}^{-1}={\left({\frac{{\mathrm{\partial }}^2{\mathrm{H}}_{\mathrm{sat}}^{\mathrm{PG}}}{\mathrm{\partial }{\varphi }_{\mathrm{PG}}^2}|}_{{\varphi }_{\mathrm{PG}}={\varphi }_{\mathrm{ss}}}\right)}^{-1},$$(18)with

$$\begin{array}{l}\frac{{\mathrm{\partial }}^2{\mathrm{H}}_{\mathrm{sat}}^{\mathrm{PG}}}{\mathrm{\partial }{\varphi }_{\mathrm{PG}}^2}\left({\varphi }_{\mathrm{PG}}\right)=P_{\mathrm{lin}}^{\mathrm{PG}}\left(1+\frac{2\pi P_{\mathrm{sat}}^{\mathrm{PG}}}{(\pi -4){\varphi }_{\mathrm{sat}}}\left(1-{\cos }^{-2}\left(\frac{\pi {\varphi }_{\mathrm{PG}}}{2{\varphi }_{\mathrm{sat}}}\right)\right)\right),\end{array}$$so that $P_{\mathrm{lin}}^{\mathrm{PG}}$ can be tuned according to

$$\begin{array}{l}{P}_{\mathrm{lin}}^{\mathrm{PG}}=\frac{(\pi -4){\varphi }_{\mathrm{sat}}}{C_{\mathrm{ss}}\left(2\pi P_{\mathrm{sat}}^{\mathrm{PG}}\left(1-{\cos }^{-2}\left(\frac{\pi {\varphi }_{\mathrm{ss}}}{2{\varphi }_{\mathrm{sat}}}\right)\right)+(\pi -4){\varphi }_{\mathrm{sat}}\right)},\end{array}$$(19)with $C_{\mathrm{ss}}=\frac{L_{\mathrm{P}}}{n_{\mathrm{P}}(0{)}^2}$ and L_P = L_C−L_leak.

4.1.3 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 frequency-dependent losses. This phenomenon is well described by a linear fractional order magnetic capacity (see [3, 34, 47–50] and [48], part 5). This is out the scope of the present work and is postponed to a follow-up paper. Here, we consider a magnetic resistance R_ec (Ω⁻¹) with magnetic impedance

$${\mathcal{T}}_{\mathrm{ec}}\left(s\right)=\frac{{\psi }_{\mathrm{ec}}\left(s\right)}{s {\varphi }_{\mathrm{ec}}\left(s\right)}=R_{\mathrm{ec}}.$$(20)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 ω_P = (R_ec C_ss)⁻¹ (Hz) and ${\tau }_{\mathrm{P}}=\frac{2\pi }{{\omega }_{\mathrm{P}}}$ (s), the resulting electrical impedance ${\mathcal{T}}_{\mathrm{P}}\left(s\right)=\frac{{\upsilon }_{\mathrm{P}}}{i_{\mathrm{C}}}$ is

$${\mathcal{T}}_{\mathrm{P}}\left(s\right)=\frac{{\upsilon }_{\mathrm{P}}\left(s\right)}{i_{\mathrm{C}}\left(s\right)}=\frac{n_{\mathrm{P}}(q_{\mathrm{D}}{)}^2}{R_{\mathrm{ec}}} \frac{s}{s+{\omega }_{\mathrm{P}}}.$$(21)

4.1.4 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 ${\mathcal{T}}_{\mathrm{P}}\left(s\right)$. For a blocked coil $\left(\frac{\mathrm{d}{q}_{\mathrm{D}}}{\mathrm{d}t}=0\Rightarrow v_{\mathcal{L}}=0\right)$, the total steady-state electrical impedance ${\mathcal{T}}_{\mathrm{C}}\left(s\right)=\frac{v_{\mathrm{I}}\left(s\right)}{i_{\mathrm{C}}\left(s\right)}$ measured at the coil terminals is given by

$${\mathcal{T}}_{\mathrm{C}}\left(s\right)=R_{\mathrm{C}}\left(1+\frac{s}{{\omega }_{\mathrm{C}}} \left(1+\frac{n_{\mathrm{P}}(q_{\mathrm{D}}{)}^2}{R_{\mathrm{C}}{R}_{\mathrm{ec}}} \frac{{\omega }_{\mathrm{C}}}{s+{\omega }_{\mathrm{P}}\left({\varphi }_{\mathrm{PG}}\right)}\right)\right).$$(22)

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 ${\mathcal{T}}_{\mathrm{C}}(i \omega )\underset{\omega \to \mathrm{\infty }}{\sim }{R}_{\mathrm{C}}\left(1+\frac{i \omega }{{\omega }_{\mathrm{C}}}\right)$. The inner bracket is the contribution of the proposed magnetic circuit.

4.1.5 Position-dependent force factor

The gyrator that restores the Lorentz force $f_{\mathcal{L}}$ with corresponding back electromotive force $v_{\mathcal{L}}$ is given by (see Eq. (B3) in Sect. B.2):

$$\left(\begin{array}{l}{v}_{\mathcal{L}}\\ f_{\mathcal{L}}\end{array}\right)=\left(\begin{array}{ll}0& -B{\mathcal{l}}_{\mathrm{C}}\\ B{\mathcal{l}}_{\mathrm{C}}& 0\end{array}\right) \left(\begin{array}{l}{i}_{\mathrm{C}}\\ {\mathrm{v}}_{\mathrm{C}}\end{array}\right),$$(23)with coil velocity ${\mathrm{v}}_{\mathrm{C}}=\frac{\mathrm{d}{q}_{\mathrm{D}}}{\mathrm{d}t}$ and ${\mathcal{l}}_{\mathrm{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 ${\mathcal{l}}_{\mathrm{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 ${\mathcal{l}}_{\mathrm{C}}:q_{\mathrm{D}}\mapsto {\mathcal{l}}_{\mathrm{C}}\left(q_{\mathrm{D}}\right)$ to model the latter phenomenon:

$${\mathcal{l}}_{\mathrm{C}}\left(q_{\mathrm{D}}\right)={{\mathcal{l}}_{\mathrm{C}}}^0\frac{1+\exp (-P_{\mathcal{l}})}{1+\exp \left(P_{\mathcal{l}} \left({\left(\frac{q_{\mathrm{D}}}{Q_{\mathcal{l}}}\right)}^2-1\right)\right)},$$(24)where ${{\mathcal{l}}_{\mathrm{C}}}^0$ is the total length of the coil, $Q_{\mathcal{l}}$ describes the overhang of the coil with respect to the magnetic path (see Fig. 13a; and Sect. 3.1.2 from [1]), and $P_{\mathcal{l}}$ is a shape parameter (see Fig. 13b).

4.2 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 ψ_M). The state is $\mathbf{x}=(x_{\mathrm{leak}},{\varphi }_{\mathrm{PG}},p_{\mathrm{M}},q_{\mathrm{D}}{)}^{\mathrm{T}}$ with the state associated with leakage flux $x_{\mathrm{leak}}=N_{\mathrm{C}} {\varphi }_{\mathrm{leak}}$, and the Hamiltonian is $\mathrm{H}\left(\mathbf{x}\right)={\mathbf{x}}^{\mathrm{T}}\mathbf{Qx}+{\mathrm{H}}_{\mathrm{sat}}^{\mathrm{PG}}\left(x_2\right)$ with $\mathbf{Q}=\frac{1}{2}\mathrm{diag}\left(\frac{1}{L_{\mathrm{leak}}},0,\frac{1}{M_{\mathrm{CDA}}},K_{\mathrm{SA}}\right)$ and ${\mathrm{H}}_{\mathrm{sat}}^{\mathrm{PG}}$ given in (16). The dissipation variable is $\mathbf{w}={\left(i_{\mathrm{C}},\frac{\mathrm{d}{q}_{\mathrm{D}}}{\mathrm{d}t},{\psi }_{\mathrm{PG}}\right)}^{\mathrm{T}}$ with linear dissipation law $\mathbf{z}\left(\mathbf{w}\right)=\mathrm{diag}(R_{\mathrm{C}},R_{\mathrm{SA}},{R_{\mathrm{ec}}}^{-1}) \mathbf{w}$. According to (14), the magnetic induction in the air gap involved in the electromechanical coupling (B3) is $b_{\mathrm{G}}=\frac{{\varphi }_{\mathrm{PG}}}{S_{\mathrm{G}}}$. The length of wire effectively subjected to the induction field is ${\mathcal{l}}_{\mathrm{C}}\left(x_4\right)$ given in (24). The number of wire turns effectively surrounding the pole piece involved in the electromagnetic coupling is n_P(x₄) 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.

4.3 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₂(t = 0) = ϕ_ss.

Eddy-current losses

The effect of the characteristic time τ_ec = 2πR_ec C_ss due to eddy-currents in the pole piece is illustrated by imposing several DC input voltages v_I, here −50 V, 0 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 τ_ec. In each case, the magnetic flux in the path is driven to a new steady state ${\varphi }_{\mathrm{PG}}\ne {\varphi }_{\mathrm{ss}}$. The long-term effects and the influence of the characteristic time τ_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 blocked-impedance. 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\mathcal{l}=\frac{{\varphi }_{\mathrm{PG}}}{S_{\mathrm{G}}}{\mathcal{l}}_{\mathrm{C}}\left(q_{\mathrm{D}}\right)$ is modulated by the coil position (same as model 0) and the magnetic flux 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.

5 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.

Conflict of interest

Author declared no conflict of interests.

Acknowledgments

This work is supported by the project ANR-16-CE92-0028, entitled Interconnected Infinite-Dimensional systems for Heterogeneous Media, INFIDHEM, financed by the French National Research Agency (ANR). Further information is available at https://websites.isae-supaero.fr/infidhem/the-project. The authors acknowledge Antonin Novak, Balbine Maillou (Laboratoire d’Acoustique de l’Université du Maine, UMR-CNRS 6613, Le Mans, France) for fruitful discussions on the loudspeaker mechanisms. The authors also acknowledge Robert Piéchaud for taking the time to complete the proofread of this paper.

---

References

Appendix A: Numerical method

In this section, we present the numerical method used in this paper for the simulation of models 0, 1 and 2. It is based on the appropriate definition of a discrete gradient [51] which restores the passive-guaranteed port-Hamiltonian structure (3) in discrete time so that numerical stability is guaranteed (see [29, 52] for details, in particular for an analysis of consistency of the method which respect to the time step, which we do not recall here).

To ensure stable simulation of stable dynamical system $\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}=\mathbf{f}\left(\mathbf{x}\right)$, many numerical schemes focus on the approximation quality of the time derivative, combined with operation of the vector field f. Here, we adopt an alternative point of view, by transposing the power balance (4) in the 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\circ \mathbf{x}$. This is the case of Euler scheme, for which first order approximation of the differential applications $\mathrm{d}\mathbf{x}(t,\mathrm{d}t)=\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}\left(t\right) \mathrm{d}t$ and dH(x, dx) = ∇H(x)^Tdx on the sample grid $t=k \delta t, k\in \mathbb{Z}$ are given by

$$\delta \mathbf{x}(k,\delta t)=\mathbf{x}(k+1)-\mathbf{x}\left(k\right),$$(A1)

$$\begin{array}{cc}\delta H\left(\mathbf{x},\delta \mathbf{x}\right)& =H\left(\mathbf{x}+\delta \mathbf{x}\right)-H\left(\mathbf{x}\right)\\ & ={\nabla }_{d}H(\mathbf{x},\mathbf{x}+\delta \mathbf{x}{)}^{\mathrm{T}} \delta \mathbf{x}.\end{array}$$(A2)

For mono-variate storage components ($H\left(\mathbf{x}\right)={\sum }_{n=1}^{n_{\mathbf{x}}}\mathrm{}{H}_n\left(x_n\right)$), the solution can be built element-wise with the n-th coordinate given by

$$\left[{\nabla }_{d}H\right(\mathbf{x},\mathbf{x}+\delta \mathbf{x}){]}_n=\{\begin{array}{ll}\frac{h_n(x_n+\delta x_n)-h_n\left(x_n\right)}{\delta x_n}& \mathrm{if} \delta x_n\ne 0,\\ h_n^{\mathrm{\prime }}\left(x_n\right)& \mathrm{otherwise} .\end{array}$$(A3)

A discrete chain rule is indeed recovered

$$\frac{\delta \mathrm{E}(k,\delta t)}{\delta t}={\nabla }^{d}H\left(\mathbf{x}\right(k),\mathbf{x}(k+1){)}^{\mathrm{T}} \frac{\delta \mathbf{x}(k,\delta t)}{\delta t}$$(A4)so that the following substitution in (3)

$$\begin{array}{ccc}\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t}\left(t\right)& \to & \frac{\delta \mathbf{x}(k,\delta t)}{\delta t}\\ \nabla \mathrm{H}\left(\mathbf{x}\right)& \to & {\nabla }^{d}H\left(\mathbf{x}\right(k),\mathbf{x}(k+1\left)\right)\end{array}$$(A5)leads to

$$\begin{array}{ccc}0& =& \mathbf{a}(k{)}^{\mathrm{T}} \mathbf{J} \mathbf{a}(k)=\mathbf{a}(k{)}^{\mathrm{T}} \mathbf{b}\left(k\right)\\ & =& \underset{\frac{\delta \mathrm{E}(k,\delta t)}{\delta t}}{\underbrace{\left[{\nabla }^d{H}^{\mathrm{T}} \frac{\delta \mathbf{x}}{\delta t}\right]\left(k\right)}}+\underset{{\mathrm{P}}_{\mathrm{D}}\left(k\right)}{\underbrace{\mathbf{z}\left(\mathbf{w}\right(k){)}^{\mathrm{T}} \mathbf{w}(k)}}-\underset{{\mathrm{P}}_{\mathrm{S}}\left(k\right)}{\underbrace{\mathbf{u}(k{)}^{\mathrm{T}} \mathbf{y}(k)}}.\end{array}$$(A6)

For pH systems composed of a collection of linear energy storing components with quadratic Hamiltonian $H_n\left(x_n\right)=\frac{x_n^2}{2C_n}$, we define $\mathbf{Q}=\mathrm{diag}(C_1\cdots C_{n_{\mathbf{x}}}{)}^{-1}$ so that the discrete gradient (A3) reads

$$\begin{array}{ccc}{\nabla }^{d}H(\mathbf{x},\mathbf{x}+\delta \mathbf{x})& =& \mathbf{Q}\left(\mathbf{x}\left(k\right)+\frac{\delta \mathbf{x}\left(k\right)}{2}\right),\end{array}$$(A7)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.

Appendix B: Recalls on magnetism

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.

B.1 Magnetic energy storage

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₀(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]):

$$\mathbf{b}={\mathbf{j}}_0\left(\mathbf{h}\right)+\mathbf{j}\left(\mathbf{h}\right)\simeq \mathbf{j}\left(\mathbf{h}\right)$$where we neglect the magnetization of vacuum so that h(b) = j⁻¹(b). The magnetic induction flux ϕ is defined as the flux of the magnetic induction field through a given surface $\mathcal{S}: \varphi \left(t\right)={\iint }_{\mathcal{S}}\mathrm{}\mathbf{b}\left(t\right) \mathrm{d}S=S b\left(t\right)$, where b(t) is the magnitude of the field b(t) that we assume constant over $\mathcal{S}$ and normal to the cross section. The magnetomotive force ψ is defined as the circulation of h along a closed b-field line $\mathcal{C}$ with length ${\mathcal{l}}_{\mathcal{C}}: \psi \left(b\right(t\left)\right)={\oint }_{\mathcal{C}}\mathrm{}h\left(b\right(t\left)\right) \mathrm{d}\mathcal{l}={\mathcal{l}}_{\mathcal{C}}\cdot h\left(b\right(t\left)\right)$, where we assume h(b(t)) constant along $\mathcal{C}$.

Energy storage

The variation of magnetic energy density (locally) stored in a sample of magnetic material is $\frac{\mathrm{d}\mathcal{E}}{\mathrm{d}t}=h\left(b\right)\frac{\mathrm{d}b}{\mathrm{d}t}$ (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 ${\mathcal{l}}_{\mathcal{C}}$ and cross section S is $\frac{\mathrm{d}E}{\mathrm{d}t}=S {\mathcal{l}}_{\mathcal{C}} h\left(b\right)\frac{\mathrm{d}b}{\mathrm{d}t}$. Rewriting the preceding relation for the magnetic flux ϕ = Sb and the magnetomotive force $\psi \left(b\right)={\mathcal{l}}_{\mathcal{C}}h\left(b\right)$ yields $\frac{\mathrm{d}E}{\mathrm{d}t}=\psi \left(\frac{\varphi }{S}\right) \frac{\mathrm{d}\varphi }{\mathrm{d}t}$. Thus, we can select ϕ as the state associated with the storage of magnetic energy with E(t) = H_mag(ϕ(t)) and we identify ${\mathrm{H}}_{\mathrm{mag}}^{\mathrm{\prime }}\left(\varphi \right)=\psi \left(\frac{\varphi }{S}\right)$:

$${\mathrm{H}}_{\mathrm{mag}}\left(\varphi \right)={\mathcal{l}}_{\mathcal{C}}{\int }_0^{\varphi }\mathrm{}h\left(\frac{\xi }{S}\right)\mathrm{d}\xi$$(B1)with the total energy variation $\frac{\mathrm{d}}{\mathrm{d}t}{\mathrm{H}}_{\mathrm{mag}}=\psi \frac{\mathrm{d}\varphi }{\mathrm{d}t}$.

B.2 Electromechanical coupling

Consider several windings of a conductive wire with section S_W, total length ${\mathcal{l}}_{\mathrm{W}}$, position q_W and velocity vector v _W = v_W e _W with constant direction e _W and magnitude ${\mathrm{v}}_{\mathrm{W}}=\frac{\mathrm{d}{q}_{\mathrm{W}}}{\mathrm{d}t}$. 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_{\mathrm{W}}={\iint }_{S_{\mathrm{W}}}\mathrm{}{\rho }_q {\mathbf{v}}_q\mathrm{d}S$ for the electric charge density ρ_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 $\mathrm{d}\mathcal{l}$ is subjected to the Lorentz force ${\mathbf{df}}_{\mathcal{L}}={\rho }_q S_{\mathrm{W}} \mathrm{d}\mathcal{l} \left({\mathbf{v}}_q+{\mathbf{v}}_{\mathrm{W}}\right)\times \mathbf{b}$. This force is orthogonal to the velocity v _q + v _W so that the associated mechanical power is $\mathrm{d}{\mathrm{P}}_{\mathcal{L}}=\mathrm{d}{\mathbf{f}}_{\mathcal{L}}\cdot ({\mathbf{v}}_q+{\mathbf{v}}_{\mathrm{W}})=0$. Integrating along the wire, one gets

$$\begin{array}{ccc}{\mathrm{P}}_{\mathcal{L}}& =& {\mathrm{v}}_{\mathrm{W}}.\underset{f_{\mathcal{L}}}{\underbrace{B {\mathcal{l}}_{\mathrm{W}} i_{\mathrm{W}}}}+i_{\mathrm{W}}.\underset{-v_{\mathcal{L}}}{\underbrace{B {\mathcal{l}}_{\mathrm{W}} {\mathrm{v}}_{\mathrm{W}}}}=0\\ & & \end{array}$$(B2)defining the Lorentz force $f_{\mathcal{L}}$ and the back electromotive force (voltage) $v_{\mathcal{L}}$. Notice the transfer is reversible and conservative in the sense that the outflow of energy from the electrical domain ${\mathrm{P}}_{\mathrm{elec}}=i_{\mathrm{W}} v_{\mathcal{L}}$ equals the inflow of the mechanical domain ${\mathrm{P}}_{\mathrm{meca}}={\mathrm{v}}_{\mathrm{W}} f_{\mathcal{L}}$. This corresponds to a gyrator with ratio $B {\mathcal{l}}_{\mathrm{W}}$:

$$\left(\begin{array}{l}{v}_{\mathcal{L}}\\ f_{\mathcal{L}}\end{array}\right)=\left(\begin{array}{ll}0& -B {\mathcal{l}}_{\mathrm{W}}\\ B{\mathcal{l}}_{\mathrm{W}}& 0\end{array}\right)\cdot \left(\begin{array}{l}{i}_{\mathrm{W}}\\ {\mathrm{v}}_{\mathrm{W}}\end{array}\right),$$(B3)with

$${\left(\begin{array}{l}{i}_{\mathrm{W}}\\ {\mathrm{v}}_{\mathrm{W}}\end{array}\right)}^{\mathrm{T}}\cdot \left(\begin{array}{l}{v}_{\mathcal{L}}\\ f_{\mathcal{L}}\end{array}\right)=0,$$(B4)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) ψ (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 υ) to the variation of the magnetic flux in the wire turn $\upsilon =\frac{\mathrm{d}\varphi }{\mathrm{d}t}$; and (ii) Ampere’s theorem that relates the mmf to the current in the wire ψ = 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:

$$\left(\begin{array}{l}{v}_{\mathrm{C}}\\ {\psi }_{\mathrm{C}}\end{array}\right)=\left(\begin{array}{ll}0& N_{\mathrm{C}}\\ N_{\mathrm{C}}& 0\end{array}\right)\left(\begin{array}{l}{i}_{\mathrm{C}}\\ \frac{\mathrm{d}{\varphi }_{\mathrm{C}}}{\mathrm{d}t}\end{array}\right).$$(B5)

Denoting by $s\in \mathbb{C}$ the Laplace variable, the correspondence between an impedance seen in the electrical domain ${\mathcal{Z}}_{\mathrm{elec}}\left(s\right)=\frac{{\upsilon }_{\mathrm{C}}\left(s\right)}{i_{\mathrm{C}}\left(s\right)}$ and its counterpart in the magnetic domain ${\mathcal{Z}}_{\mathrm{mag}}\left(s\right)=\frac{{\psi }_{\mathrm{C}}\left(s\right)}{s {\varphi }_{\mathrm{C}}\left(s\right)}=\frac{N_{\mathrm{C}}^2 i_{\mathrm{C}}\left(s\right)}{{\upsilon }_{\mathrm{C}}\left(s\right)}$ is given by

$${\mathcal{Z}}_{\mathrm{mag}}\left(s\right)=\frac{N_{\mathrm{C}}^2}{{\mathcal{Z}}_{\mathrm{elec}}\left(s\right)}.$$(B6)

As a result, if the path (P) is modeled by a magnetic capacity C_P, e.g. from the linearization of ${\mathrm{H}}_{\mathrm{mag}}^{\mathrm{\prime }}\left(\varphi \right)$ in (B1), the equivalent electrical inductance is $L_{\mathrm{C}}=N_{\mathrm{C}}^2 C_{\mathrm{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 ${\mathrm{P}}_{\mathrm{mag}}=\frac{\mathrm{d}{\varphi }_{\mathrm{C}}}{\mathrm{d}t} {\psi }_{\mathrm{C}}$ the power incoming the magnetic domain.

Appendix C: Storage functions

In this section, we precise the requirements on storage functions and detail the saturating storage functions we use in the loudspeakers models.

Requirements for storage functions

As stated in Section 2.2, the requirement for elementary storage function (Hamiltonian) $H:\mathbb{R}\supset \mathcal{X}\to {\mathbb{R}}_{+}$ are as follows:

1. It is positive semidefinite: H(x) ≥ 0, $\forall x\in \mathcal{X}$;

2. It is zero at origin: H(0) = 0;

3. It is radially unbounded: ${\lim }_{x\to \mathrm{\partial }\mathcal{X}}H\left(x\right)=+\mathrm{\infty }$ with $\mathrm{\partial }\mathcal{X}$ the boundary of domain $\mathcal{X}$ (could be {±∞};

4. Its second derivative is positive: $\frac{{\mathrm{d}}^{2}H}{\mathrm{d}{x}^2}\left(x\right)>0$, $\forall x\in \mathcal{X}$.

Requirements (1) and (2) ensure the associated energy $E=H\circ x:\mathbb{R}\ni t\mapsto H\left(x\right(t\left)\right)=E\left(t\right)$ is always positive or null; requirement (3) ensures the origin x = 0 is an attractor for the associated dynamical system from LaSalle invariance theroem; finally, requirement (4) ensures the derivative is monotonically increasing so that it invertible on its domain $\mathcal{X}$. A simple procedure to build such storage function is to integrate twice a non-vanishing, positive definite function $f:\mathcal{X}\to {\mathbb{R}}_{+}$, ${\lim }_{x\to \mathrm{\partial }\mathcal{X}}f\left(x\right)>0$, choosing the two integration constants $\{C_i{\}}_{i=\mathrm{1,2}}$ so that $c\left(x\right)={\int }_0^x\mathrm{}f\left(\xi \right)\mathrm{d}\xi +C_1$ and $H\left(x\right)={\int }_0^x\mathrm{}c\left(\xi \right)\mathrm{d}\xi +C_2$ are both null at x = 0. Note that it is not required the storage function to be symmetric with possibly $H\left(x\right)\ne H(-x)$. Some examples are given below (see also Fig. C1).

Examples of storage functions

• Quadratic: The simplest storage functions are the quadratic functions $H\left(x\right)=p \frac{x^2}{2}$, with (linear) derivative H′(x) = px.

• Polynomial: It is possible to construct polynomial storage function with only even order to ensure the requirements are satisfied, e.g. $H\left(x\right)=p_{\mathrm{lin}}\frac{x^2}{2}\left(1+p_{\mathrm{nl}}\frac{x^2}{2}\right)$, with derivative $H^{\mathrm{\prime }}\left(x\right)=p_{\mathrm{lin}}(x+p_{\mathrm{nl}}{x}^3)$.

Non-symmetric

A non symmetric storage functions can be constructed from two exponential functions as follows: $H\left(x\right)=p_{\mathrm{lin}}\left(\frac{\exp (p_{+} x)}{p_{+}}+\frac{\exp (p_{-} x)}{p_{-}}-\frac{p_{-}+p_{+}}{p_{-} p_{+}}\right)$ with derivative $H^{\mathrm{\prime }}\left(x\right)=p_{\mathrm{lin}}\left(\exp (p_{+} x)+\exp (p_{-} x)\right)$, where possibly $p_{+}\ne p_{-}$.

State saturating storage functions

In this work, the saturation effect of the suspension and the ferromagnetic path are described by the same idealized (symmetric) saturation curve c(x). It is built as the linear combination of basis functions c_lin(x) (linear behavior around the origin) and c_sat(x) (saturation effect):

$$c\left(x\right)=P_{\mathrm{lin}}\left(c_{\mathrm{lin}}\right(x)+P_{\mathrm{sat}}{c}_{\mathrm{sat}}\left(x)\right),$$(C1)

$$c_{\mathrm{lin}}\left(x\right)=x,$$(C2)

$$c_{\mathrm{sat}}\left(x\right)=\frac{4}{4-\pi } \left(\tan \left(\frac{\pi \cdot x}{2 x_{\mathrm{sat}}}\right)-\frac{\pi \cdot x}{2x_{\mathrm{sat}}}\right)$$(C3)with $c_{\mathrm{sat}}\left(x\right)\underrightarrow{x\to \pm x_{\mathrm{sat}}}\pm \mathrm{\infty }$, $\frac{\mathrm{\partial }{c}_{\mathrm{sat}}}{\mathrm{\partial }x}\left(0\right)=0$ so that c_sat(x) does not contribute around origin, and $c_{\mathrm{sat}}\left(\frac{1}{2}\right)=1$.

The corresponding Hamiltonian is obtained from

$$\mathrm{H}\left(x\right)={\int }_0^x\mathrm{}c\left(\xi \right) \mathrm{d}\xi =P_{\mathrm{lin}}\left(H_{\mathrm{lin}}\left(x\right)+P_{\mathrm{sat}}{H}_{\mathrm{sat}}\left(x\right)\right)$$(C4)with

$$\begin{array}{ccc}{H}_{\mathrm{lin}}\left(x\right)& =& \frac{x^2}{2},\\ H_{\mathrm{sat}}\left(x\right)& =& -\frac{8x_{\mathrm{sat}}}{\pi (4-\pi )} \left(\ln \left|\cos \left(\frac{\pi x}{2 x_{\mathrm{sat}}}\right)\right|+\frac{1}{2}{\left(\frac{\pi x}{2x_{\mathrm{sat}}}\right)}^2\right).\end{array}$$(C5)

This nonlinear saturating storage function proves positive definite providing the parameters (P_lin, P_sat) are positive, so that it can be used in structure (3), still preserving passivity.

Appendix D: Physical and technological parameters

Acronym d.u. stands for dimensionless unit

All Tables

All Figures