# Modeling the Effect of Eddy Currents on the Inductance of Loudspeaker Motors Using the Method of Lines

Oliver Munroe, Antonin Novak, Laurent Simon, Daniel Massicotte

► To cite this version:

Oliver Munroe, Antonin Novak, Laurent Simon, Daniel Massicotte. Modeling the Effect of Eddy Currents on the Inductance of Loudspeaker Motors Using the Method of Lines. Journal of the Audio Engineering Society, 2025, 73 (1/2), pp.63-74. 10.17743/jaes.2022.0188. hal-05395453

HAL Id: hal-05395453

https://hal.science/hal-05395453v1

Submitted on 3 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Modeling The Effect Of Eddy Currents On The Inductance Of Loudspeaker Motors Using The Method Of Lines

Oliver Munroe<fnref n="1" />, Antonin Novak², Laurent Simon²*, Daniel Massicotte<fnref n="1" />

*laurent.simon@univ-lemans.fr

## Abstract

This paper presents a way of modeling the effect of eddy currents on the inductance of a loudspeaker coil over a wide frequency bandwidth without the use of fractional derivatives or Foster/Cauer networks. The model presented is a physics-based model, which is derived from Maxwell's equations and discretized in space. The model is compared with Finite Element Analysis (FEA) simulations, and measurement data. The results show an excellent agreement with the reference data over the audio bandwidth, and for different coil positions. Using eight parameters, of which four require fitting and only two are dependent on the coil position, the model is simple to fit to reference data, and may be vectorized for real-time applications.

The archived file is not the final published version of the article O. Munroe, A. Novak, L. Simon & D. Massicotte (2025), "Modeling the Effect of Eddy Currents on the Inductance of Loudspeaker Motors Using the Method of Lines", Journal of the Audio Engineering Society. Vol. 73(1), pp. 63-74.

The definitive publisher-authenticated version is available online at https://doi.org/10.17743/jaes.2022.0188, Readers must contact the publisher for reprint or permission to use the material in any form.

<footnote n="1">¹ Laboratory of Signal and System Integration, Department of Electrical and Computer Engineering, Université du Québec à Trois-Rivières, Trois-Rivières, QC G9A 5H7, Canada
² Laboratoire d'Acoustique de l'Université du Mans (LAUM), UMR 6613, Institut d'Acoustique - Graduate School (IA-GS), CNRS, Le Mans Université, France</footnote>
# 1. INTRODUCTION

A loudspeaker is a type of voice coil actuator, where a voice coil is situated inside a magnetic gap. One or several permanent magnets provide a static magnetic field which is guided to and concentrated inside the magnetic gap by a geometry of low-carbon steel. The current flowing through the coil interacts with the static magnetic field and, via the Lorentz force, moves the coil.

The current flowing through the coil creates a time varying magnetic field which interacts with the surrounding steel, inducing eddy currents. These eddy currents generate their own magnetic field which is then coupled back into the voice coil. The consequence is that the inductance of the coil is no longer constant with frequency, but frequency dependent. The inductance has a tendency to decrease above a certain frequency, while the ohmic losses increase [1].

Various models have been published approximating the behavior of voice coil inductance taking into account the eddy currents. They can be separated into two groups:

* Group 1: Those that use only integer inductance elements.

* Group 2: Those that use fractional inductance elements.

The models in Group 1 take the form of Foster networks of pure resistors and inductors [2]. An example would be the RL₂ model used in the Klippel LSI module and published in [3]. The models in Group 2 take advantage of fractional calculus to approximate the eddy current behavior via fractional derivatives and inductors [1, 4, 5, 6]. A more detailed comparison of the two groups can be found in [7].

The model presented in this work does not belong to either of the two groups. Instead, several simplifications of Maxwell's equations for an axisymmetric geometry are used to develop a two dimensional model of the coil and steel. The resulting model is physics based, requires a small amount of physical parameters, and does not have the memory and computational requirements of fractional derivatives.

Section 1 details the physics and the various simplifications. Sections 2 & 3 detail the method used to discretize the physical equations, and the material used to evaluate the model. The results are introduced in section 4, and a discussion of the follows in section 5. Finally section 6 concludes on the work presented in this article.

# 2. State of the art

## 2.1 Loudspeaker Motor

A loudspeaker coil is coupled to a magnetic circuit comprised of low-carbon steel sections and several air gaps. Due to its symmetrical design, the loudspeaker coil and motor can be readily simplified into a two-dimensional axisymmetric problem as shown in Figure 1.

This geometry is then simplified into a magnetic circuit as shown in Figure 2. The total flux through the coil $\phi(t)$ due to the generated Magneto Motive Force $\mathcal{F}$ is split into flux kept inside the motor structure $\phi_M(t)$, and flux leaking outside the motor $\phi_L(t)$. The leakage flux path is represented by the dashed line.

For the circuit in Figure 2, $\mathcal{R}_{pole}$ is the reluctance of the pole, $\mathcal{R}_{gap}$ is the reluctance of the air gap, $\mathcal{R}_{TP}$ is the reluctance of the top plate, $\mathcal{R}_{mag}$ is the reluctance of the magnet, $\mathcal{R}_{BP}$ is the reluctance of the bottom plate, and $\mathcal{R}_{leak}$ is the reluctance of the leakage flux $\phi_L(t)$. The reluctance value for each portion may be calculated with

$$ \mathcal{R} = \frac{l}{\mu A}, \qquad (1) $$

where $l$ is the path length of the magnetic flux inside the portion, $\mu$ is the magnetic permeability of the portion, and $A$ is the surface area through which flows the flux. Converting a reluctance to an inductance can be done with

$$ L = \frac{n_{\text{Turns}}^2}{\mathcal{R}}, \qquad (2) $$
<figure><img src="image_3.png" /><figcaption>Figure 1. Loudspeaker geometry.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2. Equivalent magnetic circuit.</figcaption></figure>

where $L$ is the electrical inductance, $n_{Turns}$ is the number of coil windings, and $R$ is the magnetic reluctance. In the case of several susceptances in series, the magnetic flux $\phi$ is same for each one due to Equation 4. As the voltage developed across an inductor is proportional to $\frac{d\phi}{dt}$, this means that a circuit of series susceptances is analogous to a circuit of parallel inductances, with each inductance value calculated via Equation 2. Therefore, the magnetic reluctance circuit shown in Figure 2 could be represented by five parallel inductances $L_{pole}$, $L_{gap}$, $L_{TP}$, $L_{Mag}$ and $L_{BP}$, in series with the leakage inductance $L_{leak}$.

## 2.2 Maxwell's equations

Maxwell's equations describe the behavior of electromagnetic waves in time and space. When applied to matter, the equations may be written in terms of free current density $J_f$, electric field $E$, displacement field $D$, magnetic flux density $B$, and magnetic field strength $H$. Equations 3 to 6 illustrate Maxwell's equations applied to a macroscopic medium [8]:

$$ \nabla D = \rho_f, \qquad (3) $$

$$ \nabla B = 0, \qquad (4) $$

$$ \nabla \times E = - \frac{\partial B}{\partial t}, \qquad (5) $$
$$ \nabla \times H = J_f + \frac{\partial D}{\partial t}, \qquad (6) $$

where $\rho_f$ is the free electric charge density.

Equation 7 illustrates the relationship between the magnetic flux density $B$ and the magnetic field strength $H$ for a linear material:

$$ B = \mu H, \qquad (7) $$

where $\mu$ is the total permeability, defined as $\mu_0\mu_r$, where $\mu_0$ is the magnetic permeability of free space, and $\mu_r$ is the relative permeability of the medium. Finally for this work, the electrical field $E$ and the free current density $J_f$ are linked by the scalar electrical conductivity $\sigma$ via Ohm's law in the form [1, 9, 10],

$$ J_f = \sigma E. \qquad (8) $$

Equations 3 to 6 in the cylindrical coordinate system $(r, \theta, z)$ for the steel parts where $\rho_f = 0$ are then simplified with several assumptions. For audio frequencies the $\frac{\partial D}{\partial t}$ term in Equation 6 can be neglected [1], and we assume that due to the symmetry of the problem, $\frac{\partial}{\partial \theta} = 0$. This gives the result similar to that in [1, 10],

$$ \frac{\partial B_z}{\partial t} = \frac{1}{\sigma} \left[ \frac{\partial^2 H_z}{\partial r^2} + \frac{1}{r} \frac{\partial H_z}{\partial r} \right], \qquad (9) $$

which is the form of a diffusion equation, along the $r$ axis. In the case of a non vented pole, the $\frac{1}{r}$ in Equation 9 diverges at the point $r = 0$. As such, using L'Hopitals rule gives

$$ \lim_{r \to 0} \frac{1}{r} \frac{\partial H_z}{\partial r} = \lim_{r \to 0} \frac{\partial^2 H_z}{\partial r^2} \qquad (10) $$

thus at $r = 0$, Equation 9 becomes

$$ \frac{\partial B_z}{\partial t} = \frac{2}{\sigma} \left[ \frac{\partial^2 H_z}{\partial r^2} \right]. \qquad (11) $$

## 2.3 Relative Permeability

There are several types of relative permeability [11], depending on the state of the magnetic material. For this work, we will only consider the differential relative permeability, defined as

$$ \mu_{r\Delta} = \frac{1}{\mu_0} \frac{dB}{dH}. \qquad (12) $$

The differential permeability is used when there is a static magnetic field creating an operating point inside the magnetic material, superimposed with a smaller time varying excitation [11]. This matches the configuration of a loudspeaker motor, where the magnet provides the static magnetic operating point and the time varying field generated by the coil works around this point. Figure 3 shows $\mu_{r\Delta}$ referenced to the anhysteretic $B(H)$ curve for an AISI 1010 steel, with the data tips at the average motor operating point inside the pole.

# 3. Methods

## 3.1 Inductance model

For the model presented in this work, the magnetic circuit in Figure 2 is simplified again to two series resistances, $R_{mag}$ and $R_{pole}$, in parallel with the $R_{leak}$ reluctance. The structure of the electrical model of the voice coil is shown in Figure 4, where $L_{mag}$ is the inductance associated with $R_{mag}$, the $u_{pole}$ voltage source represents the contribution from eddy currents in the motor, which are lumped together with the ($R_{pole}$) reluctance, $R_e$ is the DC resistance of the coil, and $L_{leak}$ is the inductance due to $R_{leak}$. $L_{mag}$ and $L_{leak}$ are considered independent of
<figure><img src="image_5.png" /><figcaption>**Figure 3.** AISI 1010 Steel B(H) curve with $\mu_{r\Delta}$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Voice coil model.</figcaption></figure>

both frequency and coil current. As the coil moves out of the motor the coupling between the coil and the motor's magnetic circuit will decrease resulting in $n_{Turns}$ and $L_{mag}$ being dependent on the coil position $x$, which is in the $z$ direction. It should be noted here that the $n_{Turns}(x)$ parameter is bounded, $0 \le n_{Turns}(x) \le c_{Turns}$, where $c_{Turns}$ is the total number of coil turns and is a constant.

The model structure bears a resemblance to the model in [4], and more closely of that presented in [6], the main difference being the way the influence of the pole is approximated. In [4, 6], the influence of the eddy currents in the pole is represented by a fractional inductance, whereas for this work the influence of the eddy currents in the pole is represented by a discretization in space of Equation 9 and considers the relative permeability to be constant for all values of $H_z$.

### 3.1.1 Assumed spatial dependence of $n_{Turns}(x)$ & $L_{Mag}(x)$

The $n_{Turns}$ and $L_{Mag}$ parameters may have different spatial variations due to the geometry of the motor. For the $n_{Turns}$ parameter, a proposed approximation is to observe the amount of coil above the pole and subtract the equivalent amount of turns. For example at the coil's rest position (0 mm), the top of the coil is approximately 0.7 mm above the top of the pole. As the coil has two layers of wire of radius $r_w = 0.15$ mm, this is equivalent to around 5 turns, thus the theoretical $n_{Turns}$ is $78 - 5 = 73$. The same logic is applied for the $L_{Mag}$ parameter,
except it is assumed that only the turns below point A in Figure 1 are coupled to the magnetic circuit of the motor. Thus the theoretical values are generated using Equations 1 and 2 to compute $L_{Mag}$ using the amount of turns under the point A of the top plate. It should be noted here that this is an oversimplification of the physics, and a much more in-depth analysis can be found in [12].

## 3.2 Pole model

The influence of the eddy currents inside the steel pole is represented by the $u_{pole}$ voltage source, which is the result of the partial differential equation in Equation 9, where the time derivative of $B_z$ depends on the spatial gradient of $H_z$. Other works [1, 4, 5, 6, 10, 13] use a fractional derivative operator in order to approximate this relation. However in this work the Method Of Lines [14], a method for solving partial differential equations, is used to solve Equation 9.

To this end, the pole is discretized into $N$ points along the $r$ axis (with $N \gg 1$), where the first point ($n = 0$) is on the outer radius of the pole and last point ($n = N - 1$) is on the axis of symmetry. In this case, the spatial resolution $\Delta r$ is defined as

$$ \Delta r = \frac{r_{pole}}{N - 1}, \qquad (13) $$

and the discretized points as

$$ r(n) = r_{pole} - n\Delta r, \qquad (14) $$

for $n \in [0, N - 1]$.

For each point $n \in [1, N - 1]$ in space along $r$ we may use the approximations

$$ \frac{\partial H_z^n}{\partial r} \approx \frac{H_z^n - H_z^{n-1}}{\Delta r} := \frac{\Delta H_z^n}{\Delta r} \qquad (15) $$

and

$$ \frac{\partial^2 H_z^n}{\partial r^2} \approx \frac{H_z^{n+1} - 2H_z^n + H_z^{n-1}}{\Delta r^2} := \frac{\Delta^2 H_z^n}{\Delta r^2}. \qquad (16) $$

The end result is a system of $N - 1$ ordinary differential equations of the form

$$ \frac{\partial B_z^n}{\partial t} \approx \frac{1}{\sigma \Delta r^2} \left[ \Delta^2 H_z^n + \Delta r \frac{\Delta H_z^n}{r(n)} \right]. \qquad (17) $$

The quantity $\frac{\partial B_z^n}{\partial t}$ is then integrated using a numerical integration scheme giving $B_z$, which is then divided by $\mu = \mu_0\mu_{r\Delta}$ to give $H_z$. The input is modeled as the magnetic field strength at the surface of the material $H_z^0$, and is a function of the input current $i_1(t)$ and a gain factor $K_1$, described later in this work, $H_z^0 = K_1 i_1(t)$. Equations 18 to 22 give a short example how the time stepping works when the initial conditions for the $H_z$ and $B_z$ vectors are 0, as used in this work. At the first time step $t = \Delta t$, we have

$$ \frac{\partial B_z^1}{\partial t}(\Delta t) \approx \frac{1}{\sigma \Delta r^2} \left[ \Delta^2 H_z^1(\Delta t) + \Delta r \frac{\Delta H_z^1(\Delta t)}{r(1)} \right], \qquad (18) $$

with

$$ \Delta^2 H_z^1(\Delta t) = H_z^2(\Delta t) - 2H_z^1(\Delta t) + H_z^0(\Delta t), \qquad (19) $$

and

$$ \Delta H_z^1(\Delta t)) = H_z^1(\Delta t) - H_z^0(\Delta t). \qquad (20) $$

Then, using the Euler forward method

$$ B_z^1(2\Delta t) = B_z^1(\Delta t) + T_s * \frac{\partial B_z^1}{\partial t}(\Delta t), \qquad (21) $$
where $T_s$ is the sample rate. Finally,

$$H_z^1(2\Delta t) = \frac{B_z^1(2\Delta t)}{\mu}. \qquad (22)$$

This process repeats every time step and the values of $B_z^n$ and $H_z^n$ propagate through the vectors depending on the physics.

The voltage $u_{pole}(t)$ generated due to the influence of the pole can be calculated by integrating $\frac{\partial B_z(t)}{\partial t}$ across the surface of the pole to give the total time varying magnetic flux inside the coil, thus for a single turn:

$$u_{pole}(t) = \int_{0}^{2\pi} \int_{0}^{r_{pole}} \frac{\partial B_z(t)}{\partial t} r dr d\theta, \qquad (23)$$

which may also be written

$$u_{pole}(t) \approx 2\pi\Delta r \sum_{n=0}^{N-1} \frac{\partial B_z^n(t)}{\partial t} r(n). \qquad (24)$$

### 3.2.1 Axis of symmetry

The point at $n = N - 1$ is situated on the axis of symmetry, therefore the points on either side of this boundary will be equal, giving $H_z^N = H_z^{N-2}$, where $H_z^N$ is a fictitious quantity arising from the symmetry of the problem. By using this result with Equations 11 and 16, we arrive at the expression for the point situated on the axis of symmetry at $n = N - 1$

$$\frac{\partial B_z^{N-1}}{\partial t} = \frac{4}{\sigma} \left[ \frac{H_z^{N-2} - H_z^{N-1}}{\Delta r^2} \right]. \qquad (25)$$

### 3.2.2 Vented Pole

In the case of a vented pole, or where the pole geometry does not extend to the axis of symmetry, we apply the Neumann boundary condition [15] at the pole/air boundary

$$\frac{\partial H_z^{N-1}}{\partial r} = 0, \qquad (26)$$

which, using the spatial discretization in Equation 15 gives

$$H_z^{N-1} = H_z^{N-2}. \qquad (27)$$

It should be noted here that this boundary condition is an approximation inspired by physics and found to work well when used with the pole model. In reality, the value of $H_z$ is discontinuous across the boundary and dependent on the surface current at the boundary [16]. However, this error is only visible at very low frequencies where eddy currents flow through the whole pole. As the frequency of $H_z$ increases, the skin effect pushes $H_z$ towards the pole's outer radius, reducing the effect of the boundary at the pole's inner radius.

### 3.2.3 Skin depth and spatial resolution

The voltage $u_{pole}$ depends on the spatial gradient of the magnetic field strength $H_z$. Therefore, the $u_{pole}$ will drop to zero above a certain frequency where the skin depth is smaller than the spatial resolution of the pole model, shorting out the $L_{Mag}$ inductance. Of course it is possible to increase the number of points, however the system of equations may become unstable if the sampling frequency is too low. For example when using the explicit one-step Euler Forward method, if $\Delta r$ is decreased by a factor of two then the minimum sampling frequency to guarantee stability is increased by a factor of four due to the $\frac{1}{\Delta r^2}$ factor [14].
# 4. Material

In order to perform an in-depth analysis of the model proposed in this work, both FEA and measurement data is used. A large amount of data and hardware was already available due to work done in [7], which is partly reused here. The analysis is limited to the frequency dependency due to eddy currents, and the spatial variation of the parameters due to the coil moving inside the motor. Analysis and modeling of the distortion due to the reluctance force, the B(H) curve, and the magnetic hysteresis are outside the scope of this article.

## 4.1 Drive units

The drive unit used to assess the model accuracy is the Faital Pro 5FE100. A large amount of measurement and simulation data was already available due to the work done in [7] and [17].

## 4.2 Measured Data

The measurements are performed using a specially designed jig [17], where the voice coil is blocked at a specific position x. For the electrical impedance measurements, the voice coil is then stimulated with a 60 second exponential sweep using a voltage source, with the sweep parameters described in Table 1. The 0.1 Hz start frequency is due to the window function applied to the signal, and the useful sweep range starts at 1 Hz. The measurements are repeated for positions between -8 mm and + 8 mm with a step of 2 mm. The DC resistance $R_e$ is measured separately before each measurement to account for any changes due to thermal effects.

The voice coil voltage and current are measured using hardware detailed in [7], comprised of INA128 instrumentation operational amplifiers. The voice coil voltage and current signals are acquired using a Data Translation DT9837C Data Acquisition System, and post processed in Python and Matlab. The frequency domain data from the measurements is cropped to below 1 kHz as the resonant behavior of the jig is visible in the results above this frequency.

<table><caption>**Table 1.** Measurement Sweep Parameters.</caption><thead><tr><td>Parameter</td><td>Value</td><td>Unit</td><td>Description</td></tr></thead><tbody><tr><td>F<sub>s</sub></td><td>96</td><td>kHz</td><td>Sampling frequency</td></tr><tr><td>f<sub>1</sub></td><td>0.1</td><td>Hz</td><td>Sweep start frequency</td></tr><tr><td>f<sub>2</sub></td><td>48</td><td>kHz</td><td>Sweep end frequency</td></tr><tr><td>T</td><td>60</td><td>s</td><td>Sweep duration</td></tr><tr><td>A</td><td>3</td><td>V</td><td>Sweep amplitude</td></tr></tbody></table>

## 4.3 Simulated Data

The FEMM software [18] is used to generate the simulated impedance data from 2 Hz up to 20 kHz, for the coil positions between -10 mm and +10 mm, with a 2 mm step. For each frequency point in the study, the FEA data is generated using an AC analysis with frozen permeability based on an initial static solution. The geometry of the model is generated by measuring the dimensions of the actual device, and the material parameters chosen to best fit the measured electrical impedance at a selection of coil positions.

## 4.4 Model Parameters

The material parameters used in the model are shown in Table 2, and the physical dimensions of the various motor components used in this work are shown in Table 3. A static FEA simulation of the motor shows that the average $B_z$ value in the pole due to the permanent magnet is 1.4 T. Using the data in Figure 3 provides us with a value for $\mu_{r\Delta}$.
<table><caption>**Table 2.** Material Parameters.</caption><thead><tr><td>Parameter</td><td>Value</td><td>Unit</td></tr></thead><tbody><tr><td>σ</td><td>5.8</td><td>MS·m<sup>-1</sup></td></tr><tr><td>μ<sub>0</sub></td><td>4π × 10<sup>-7</sup></td><td>N·A<sup>-2</sup></td></tr><tr><td>μ<sub>rΔ</sub></td><td>300.9</td><td>-</td></tr></tbody></table>

<table><caption>**Table 3.** Dimensions.</caption><thead><tr><td>Parameter</td><td>Value</td><td>Unit</td><td>Description</td></tr></thead><tbody><tr><td>c<sub>Turns</sub></td><td>78</td><td>-</td><td>Coil turns</td></tr><tr><td>l<sub>coil</sub></td><td>12</td><td>mm</td><td>Coil length</td></tr><tr><td>r<sub>w</sub></td><td>0.15</td><td>mm</td><td>Wire radius</td></tr><tr><td>r<sub>pole</sub></td><td>15.9</td><td>mm</td><td>Pole outer radius</td></tr><tr><td>r<sub>coil</sub></td><td>16.3</td><td>mm</td><td>Coil radius</td></tr><tr><td>r<sub>o<sub>mag</sub></sub></td><td>45</td><td>mm</td><td>Magnet outer radius</td></tr><tr><td>r<sub>i<sub>mag</sub></sub></td><td>22</td><td>mm</td><td>Magnet inner radius</td></tr><tr><td>h<sub>mag</sub></td><td>18</td><td>mm</td><td>Magnet thickness</td></tr></tbody></table>

The physical parameters in Tables 2 and 3 are then used to define the model parameters, described in Table 4. Here, N is chosen as a power of two for illustrative purposes only, one may use an arbitrary number of points.

<table><caption>**Table 4.** Model Parameters.</caption><thead><tr><td>Parameter</td><td>Value</td><td>Description</td></tr></thead><tbody><tr><td>N</td><td>128</td><td>Number of spatial points</td></tr><tr><td>Δr</td><td>125 μm</td><td>Spatial resolution</td></tr><tr><td>L<sub>mag</sub></td><td>2.05 mH</td><td>Inductance of Magnet</td></tr><tr><td>L<sub>leak</sub></td><td>0.24 mH</td><td>Leakage inductance</td></tr><tr><td>K<sub>1</sub></td><td>9.6 m<sup>-1</sup></td><td>Conversion from current i<sub>1</sub>(t) to H<sub>0</sub></td></tr><tr><td>K<sub>2</sub></td><td>10.99 Ω·m<sup>-2</sup></td><td>1/σ(Δr<sup>2</sup>)</td></tr><tr><td>K<sub>3</sub></td><td>7.86 × 10<sup>-4</sup> m</td><td>2πΔr</td></tr></tbody></table>

The $L_{mag}$ parameter is calculated using Equations 1 and 2, with $c_{Turns}$, and $μ_r = 1$. The $L_{leak}$ parameter is calculated using Wheeler's formula for an air cored coil [19, 20],

$$L_{leak} [H] = \frac{31.6 \times c_{Turns}^2 \times r_{coil}^2}{6r_{coil} + 9l_{coil} + 40r_w} 10^{-6}. \quad (28)$$

An initial value for the $K_1$ coefficient can be found using

$$K_1 = \frac{1}{l_{coil}} \times \frac{2}{\sqrt{\mu_{r\Delta}}}. \quad (29)$$

The $\frac{1}{l_{coil}}$ factor is from Ampère's circuital law applied to a single turn, and the $\frac{2}{\sqrt{\mu_{r\Delta}}}$ factor comes from the transmission of an electromagnetic wave traveling from one medium to another. The calculations for this factor are detailed in the Appendix A.0.1.
## 4.5 Model Implementation

In this work, the blocked coil behavior of the voice coil inside the motor is represented by the following differential equations

$$ \frac{di(t)}{dt} = \frac{u_{stim}(t) - R_e i(t) - u_{pole}(i_1(t)) n_{Turns}}{L_{leak}}, \quad (30) $$

and

$$ \frac{di_2(t)}{dt} = n_{Turns} \frac{u_{pole}(i_1(t))}{L_{Mag}}, \quad (31) $$

where $i_1(t) = (i(t) - i_2(t))n_{Turns}$ and $u_{stim}(t)$ is the stimulus. The pole is modeled as extending all the way to the axis of symmetry. Pseudo code is provided in Appendix A.0.2.

## 4.6 Parameter Fitting

The fitting is done in two steps. First, the parameters $L_{leak}$, $\mu_{r\Delta}$, and $L_{Mag}$ are fitted to the simulated electrical impedance of the coil fully inside the motor (coil position at -10 mm). At this position, the coupling between the coil and the motor is assumed to be at its maximum and thus the $n_{Turns}$ parameter may be fixed at the total number of turns of the coil. Once these three parameters are fitted, then the spatial variations of the $n_{Turns}$ and $L_{Mag}$ parameters are fitted to the simulated electrical impedance of the coil at different positions in the motor.

The coil model is run in the time domain and stimulated using a swept-sine. The modeled electrical impedance $Z$ was then obtained via Fast Fourier Transforms of the modeled coil voltage and current. The quantity used to fit the model is the apparent inductance $L_{app}$ defined as

$$ L_{app} = \frac{\Im(Z)}{\omega}, \quad (32) $$

with $\omega$ the angular frequency defined as $\omega = 2\pi f$, and $\Im(Z)$ the imaginary part of the electrical impedance $Z$.

The sweep was generated using the method in [21] and with the parameters in Table 5. Again, the 0.1 Hz start frequency is due to the window function applied to the signal, and the useful sweep range starts at 1 Hz. A 10 second sweep is used instead of 60 seconds in order to speed up the fitting time, as the model is run in the time domain. Finally, the sweep end frequency $f_2$ is set at the Nyquist frequency for noise reasons.

The reference data comes from FEA simulations, using the FEMM software [18]. The use of FEA to provide the reference data is advantageous as there are no thermal effects or mechanical resonances. It was shown in [7] and [17] that the FEA simulation can be used to accurately model the linear electrical impedance of the voice coil.

<table><caption>Table 5. Sweep Parameters.</caption><thead><tr><td>Parameter</td><td>Value</td><td>Unit</td><td>Description</td></tr></thead><tbody><tr><td>F<sub>s</sub></td><td>200</td><td>kHz</td><td>Sampling frequency</td></tr><tr><td>f<sub>1</sub></td><td>0.1</td><td>Hz</td><td>Sweep start frequency</td></tr><tr><td>f<sub>2</sub></td><td>100</td><td>kHz</td><td>Sweep end frequency</td></tr><tr><td>T</td><td>10</td><td>s</td><td>Sweep duration</td></tr><tr><td>A</td><td>3</td><td>V</td><td>Sweep amplitude</td></tr></tbody></table>

# 5. Results

## 5.1 $L_{leak}$, $\mu_{r\Delta}$, and $L_{Mag}$

The initial guess for the $L_{leak}$, $\mu_{r\Delta}$, and $L_{Mag}$ parameters are compared to the result of the parameter fitting in Table 6, and the final fitted curves for $L_{app}$ and $R_{app}$, ($R_{app} = \Re(Z) - R_e$) are shown in Figure 5. The initial guess is very close to the final result, and the model shows acceptable accuracy across the audio bandwidth.
<table><caption>Table 6. Parameter Fitting.</caption><thead><tr><td>Parameter</td><td>Initial Value</td><td>Fitted Value</td></tr></thead><tbody><tr><td>L<sub>leak</sub></td><td>0.24 mH</td><td>0.36 mH</td></tr><tr><td>μ<sub>rΔ</sub></td><td>300.9</td><td>299</td></tr><tr><td>L<sub>Mag</sub></td><td>2.05 mH</td><td>2.1 mH</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>Figure 5. Comparison between the model and the reference FEA simulation.</figcaption></figure>

## 5.2 Spatial variations of $n_{Turns}$ & $L_{Mag}$

Figures 6 and 7 show the spatial variations of the $n_{Turns}$ and $L_{Mag}$ parameters resulting from fitting the model parameters to the reference data over several coil positions. The model was fitted to reference data at coil positions ranging from -10 mm to + 10 mm with a 2 mm step. As expected, the $n_{Turns}$ and $L_{Mag}$ parameters decrease as the coil moves out of the gap. Polynomial functions are used to represent the trends of the spatial variations of each parameter. The trends of the $n_{Turns}$ and $L_{Mag}$ parameters are represented by 7th and 3rd order polynomial functions respectively.

Figure 8 compares the model to the reference FEA data used for fitting and also the measurement data for five different coil positions (-8 mm, -4 mm, 0 mm, 4 mm, 8 mm). Once again the model shows very good accuracy over a wide bandwidth for each coil position. The measurement data only goes up to 1 kHz as above this frequency the mechanical modes of the measurement jig influenced the results.

## 5.3 Effect of N

Figure 9 shows how the accuracy of the model is affected by the number of spatial discretization points, showing the effect on the apparent inductance and resistance at a coil position of - 10 mm. The model parameters used are those fitted to the model with 128 points, and as such do not change for each value of N. As the number of points is reduced, the bandwidth over which the model is accurate is decreased. Moving from 128 points to 64 points effects the accuracy of the model of a wide frequency range, due to the skin depth. Using the parameter values presented
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Spatial variation of the *n*<sub>Turns</sub> parameter. Circles show the fitted values, crosses show the theoretical values.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Spatial variation of the *L*<sub>Mag</sub> parameter. Circles show the fitted values, crosses show the theoretical values.</figcaption></figure>

in this work for σ and μ, Table 7 gives the skin depth for several frequencies.

The time to run the model was found not to depend on the amount of *N* points used during this study. For each simulation for *N* ∈ [8, 128], the simulation time was 12 seconds for the 10 s stimulus running at a sampling frequency of 200 kHz.

## 6. Discussion
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Comparison between the model, a FEA simulation, and a measurement performed on a drive unit (DUT). Results shown at -8 mm, -4 mm, 0 mm, +4 mm and +8 mm.</figcaption></figure>

<table><caption>**Table 7.** Skin Depth in mm.</caption><thead><tr><th>30 Hz</th><th>130 Hz</th><th>500 Hz</th><th>3 kHz</th><th>10 kHz</th></tr></thead><tbody><tr><td>2.2</td><td>1.07</td><td>0.55</td><td>0.22</td><td>0.128</td></tr></tbody></table>

## 6.1 Model accuracy

Observation of the results presented in Figures 5 and 8 proves that the model is very accurate over most of the audio bandwidth and for each coil position. However, the number of spatial discretization points has a very large influence on the accuracy of the model, as shown in Figure 9. If the spatial resolution is lower than the skin depth, the model cannot correctly represent the spatial variations of $H_z$, and the $u_{pole}$ voltage drops to zero. Once the $L_{mag}$ inductance is shorted out by $u_{pole}$, the only inductance remaining is the $L_{leak}$ inductance. By comparing Table 7 and Figure 9, it is possible to observe how the spatial resolution of 64 points is not sufficient to resolve the eddy currents for frequencies above approximately 3 kHz. The error appears gradually, partly due to the skin effect and partly to the parallel association of $L_{Mag}$ and $u_{pole}$.

For the model parameters used in this work the skin depth $\delta$ at 20 kHz is approximately 87 $\mu$m. For acceptable accuracy over the audio bandwidth, a spatial resolution of approximately 125 $\mu$m seems to be sufficient. In theory, 183 points would be required to give a spatial resolution of 87 $\mu$m, however above 10 kHz, the apparent inductance is mostly due to the leakage inductance as the eddy currents have effectively shielded all the steel parts. As such, the error induced by the spatial resolution is less visible.

Using a larger number of points would be preferable, however as described in Section 2.2.3, doubling the number of points requires a quadrupling of the minimum sampling frequency in order to guarantee stability. For example, with the parameters in this work 256 points required increasing the sampling frequency to 300 kHz, but did not significantly improve the accuracy.

Thus the results shown here represent a compromise between accuracy and run time, it is up to the reader to decide which compromise best fits their situation.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Comparison between the model and a FEA simulation for different numbers of spatial discretization points.</figcaption></figure>

## 6.2 Model parameters

There are a total of eight parameters required for this model: $L_{Mag}$, $L_{leak}$, $n_{Turns}$, $\mu_{r\Delta}$, $\sigma$, $N$, $r_{pole}$, and $l_{Coil}$. Only four of the parameters require fitting, with the others being defined by the physical aspects of the motor and coil. Of these four, only two parameters are dependent on coil position: $L_{Mag}$ and $n_{Turns}$. Due to the physics-based modeling, initial values for all the parameters may be calculated using the physical aspects of the drive unit motor.

As shown in Table 6, the initial values are very close to the fitted values. The $L_{leak}$ parameter is believed to be higher than the initial value due to the spatial resolution of the pole model, and an oversimplification of the physics. As the skin depth decreases to values smaller than the spatial resolution of the pole model, the $u_{pole}$ voltage will tend to 0, effectively shorting out the $L_{mag}$ inductance. This will lead to an underestimation of both the apparent inductance and resistance above a certain frequency. As the fitting is done with a time domain model, these effects are taken into account when optimizing the parameters, leading to an increase in the leakage inductance to compensate for the insufficient spatial resolution. There is scope to improve the fitting procedure, and improve the model, and work will continue on both these points.

The spatial variations of both $L_{Mag}$ and $n_{Turns}$ are found to be justified from a physics point of view. The $n_{Turns}$ parameter starts at the total number of coil turns when the coil is fully immersed inside the motor system, and then decreases as the coil moves above the rest position (0 mm) and out of the motor system. We can see that the assumed values for $n_{Turns}$ provide a good initial guess. The differences may come from a number of issues, one of which being that the theoretical values consider the coupling to be binary; a turn is either coupled to the pole or not.

The $L_{Mag}$ parameter also decreases as the coil moves up and out of the motor system, but much faster than the $n_{Turns}$ parameter. As the theoretical values for $L_{Mag}$ show a sharper decrease than the fitted values, it is thought that the coil stays coupled to the magnetic circuit of the motor until the top of the coil is level with point B in Figure 1. The predicted spatial variation of the inductance value in [12] for the overhung coil shows a very similar trend to the spatial variation of the $L_{Mag}$ parameter. It should therefore be possible to use the work in [12] to simplify the fitting of the parameters for different coil positions.
Parameters $n_{Turns}$ and $L_{Mag}$ do not share the same spatial variation, this is thought to be due to the geometry of the motor, where the pole extends above the top plate. This means that as the coil moves out of the magnetic circuit of the motor, it will still be relatively well coupled to the pole. This can be clearly seen in Figure 6, where the $n_{Turns}$ parameter is relatively constant up until the 0 mm position. As the coil moves up and away from the 0 mm position, the coupling ($n_{Turns}$) decreases. In comparison, Figure 7 shows that the $L_{Mag}$ parameter decreases as soon as the coil moves up and out of the motor. For coil positions above +6 mm, the $L_{Mag}$ parameter can be seen to tend towards a constant value. This can be explained physically by the fact that once outside the influence of the motors magnetic circuit, the inductance will revert to that of an air cored coil.

## 6.3 Model simplicity

The model presented here requires eight parameters ($L_{Mag}$, $L_{leak}$, $n_{Turns}$, $\mu_{r\Delta}$, $\sigma$, $N$, $r_{pole}$, $l_{Coil}$), for which the initial values may be calculated using the physical characteristics of the motor and coil. Only two parameters depend on the coil position, and the parameter describing the $B(H)$ relationship is independent on the coil position. The pole model is a system of N-1 ordinary differential equations, where N defines the number of points in space. As can be seen in the results presented in Figures 5 and 8, $N = 128$ points is sufficient to represent the behavior of the voice coil inductance across most of the audio bandwidth.

In comparison:

* A fractional derivative using the Grünwald-Letnikov approximation usually [22] requires several thousand memory points that need to be summed at each time step, and may have weights that depend on the coil position [13] which requires a recalculation of the FIR weights using either the recurrence relation method or the FFT method [22] for different coil positions.

* A $RL_N$ network requires several cells to represent the behavior across the same bandwidth. This leads to a decrease in the $L_N$ parameter value for each cell thus increasing the stiffness of the equations, and an increase in the number of parameters to fit.

Both the fractional derivative and $RL_N$ type models either require more parameters or more calculations, and the link between the physical system and the model parameters may not be as clear as the model presented here. Finally, the structure of the pole model equations means they can be vectorized, and each ODE calculated in parallel.

Using Matlab 2023b on a windows PC with an i7 9850H processor and 16 Gb RAM, the model script took around 12 seconds to run using the 10 s stimulus at a sampling frequency of 200 KHz, for $N \in [8, 128]$. In this range, the number of points has little effect on the run time due to the vectorized equations. When using more points however, it is important to remember that doubling the number of points requiring a quadrupling of the minimum sampling frequency to guarantee stability due to the $\frac{1}{\Delta r^2}$ factor, which may increase the run time of the model.

## 7. Conclusion

The model presented in this work is a physics-based model with a small (8) amount of parameters. Two parameters are dependent on the voice coil position inside the motor. The model is based upon the application of Maxwell's equations to an axisymmetric geometry, discretized in space and solved using traditional numerical integration. Most, if not all, motor topologies can be reduced to a magnetic circuit that is identical, or very similar to the in Figure 2, and as such the model presented in this work should be applicable to most motor topologies.

As it is a physics-based model, the initial guess for parameter fitting may be calculated using the motor dimensions and material parameters.

With a spatial resolution of ≈ 125 µm, the model is capable of correctly predicting the electrical impedance of a loudspeaker coil across the audio bandwidth and for positions inside and outside of the motor.
The structure of the model means it can be parallelized for more efficient computations in real time.

## References

[1] J. Vanderkooy, "A model of loudspeaker driver impedance incorporating eddy currents in the pole structure," *J. Audio Eng. Soc.*, vol. 37, no. 3, pp. 119-128, March 1989.

[2] R. M. Foster, "A reactance theorem," *Bell System Technical Journal*, vol. 3, no. 2, pp. 259-267, April 1924.

[3] M. Dodd, W. Klippel, and J. Oclee-brown, "Voice coil impedance as a function of frequency and displacement," *J. Audio Eng. Soc.*, October 2004.

[4] I. Anazawa, "Impact of the coupling factor on lossy voice coil impedance," in *Audio Engineering Society Convention 146*, March 2019.

[5] W. M. Leach, jr., "Loudspeaker voice-coil inductance losses: Circuit models, parameter estimation, and effect on frequency response," *J. Audio Eng. Soc.*, vol. 50, no. 6, pp. 442-450, June 2002.

[6] K. Thorborg and C. Futtrup, "Electrodynamic transducer model incorporating semi-inductance and means for shorting ac magnetization," *J. Audio Eng. Soc.*, vol. 59, no. 9, pp. 612-627, September 2011.

[7] O. Munroe, "Real time loudspeaker control," Ph.D. dissertation, November 2022, le Mans 2022. [Online]. Available: www.theses.fr/2022LEMA1022

[8] J. D. Jackson, *Classical electrodynamics*, 3rd ed. John Wiley & Sons, 1999.

[9] O. Darrigol, *Electrodynamics from Ampère to Einstein*. Clarendon Press, 2000.

[10] A. Bezzola, P. Brunet, and S. Yuan, "Variable fractional order analysis of loudspeaker transducers: Theory, simulations, measurements, and synthesis," in *Audio Engineering Society Convention 143*, October 2017.

[11] B. Cullity and C. Graham, *Introduction to Magnetic Materials*. Wiley, 2009.

[12] F. T. Agerkvist and H. Franz, "On the interdependence of loudspeaker motor nonlinearities," in *Audio Engineering Society Convention 145*, October 2018.

[13] A. W. King and F. T. Agerkvist, "Position dependence of fractional derivative models for loudspeaker voice coils with lossy inductance," in *Audio Engineering Society Convention 142*, May 2017.

[14] W. E. Schiesser and G. W. Griffiths, *A Compendium of Partial Differential Equation Models: Method of Lines Analysis with Matlab*. Cambridge University Press, 2009.

[15] A. H.-D. Cheng and D. T. Cheng, "Heritage and early history of the boundary element method," *Engineering Analysis with Boundary Elements*, vol. 29, no. 3, pp. 268-302, March 2005.

[16] W. Hayt, *Engineering Electromagnetics*, ser. Electrical Engineering Series. McGraw-Hill, 1989.

[17] O. Munroe, A. Novak, and L. Simon, "Reluctance force modeling and compensation," *J. Audio Eng. Soc.*, vol. 70, no. 3, pp. 177-184, March 2022.

[18] D. C. Meeker, "Finite element method magnetics." [Online]. Available: https://www.femm.info

[19] H. A. Wheeler, "Simple inductance formulas for radio coils," *Proceedings of the Institute of Radio Engineers*, vol. 16, pp. 1398-1400, October 1928.
[20] —, “Formulas for the skin effect,” Proceedings of the IRE, vol. 30, pp. 412–424, September 1942.

[21] A. Novak, P. Lotton, and L. Simon, “Synchronized swept-sine: Theory, application, and implementation,” J. Audio Eng. Soc., vol. 63, no. 10, pp. 786–798, October 2015.

[22] P. Igor, *Fractional Differential Equations*. Academic Press, 1999.

# APPENDIX

## A.0.1 Transmission into steel

In air there are no currents ($J = 0$), and the charge density is low. Maxwell's equations can thus be reduced as shown in Equations A.1 to A.4

$$ \nabla E = 0, \qquad (A.1) $$

$$ \nabla B = 0, \qquad (A.2) $$

$$ \nabla \times E = -\frac{\partial B}{\partial t}, \qquad (A.3) $$

$$ \nabla \times B = \mu\epsilon\frac{\partial E}{\partial t}. \qquad (A.4) $$

Combining Equations A.3 and A.4 gives

$$ \nabla^2 B = \mu\epsilon\frac{\partial^2 B}{\partial t^2}. \qquad (A.5) $$

Equation A.5 has the form of a standard wave equation. Thus the magnetic field emitted by the coil will obey the physics of a wave. When the wave arrives at the steel pole, some of the energy will be reflected back, while some will be transmitted. The incident wave is assumed to be planar and normal to the surface of the pole. Thus, using Fresnel's equations, it is possible to find the reflection coefficient, $K_R$, describing the amount of the incident wave reflected back to the coil. For a normal incident wave, the reflection coefficient is

$$ K_R = \frac{n_1 - n_2}{n_1 + n_2}, \qquad (A.6) $$

with $n_1 = \sqrt{\frac{\mu_1}{\epsilon_1}}$, $n_2 = \sqrt{\frac{\mu_2}{\epsilon_2}}$. The transmission coefficient $K_T$ is defined as $1 - K_R$, giving

$$ K_T = \frac{2n_1}{n_1 + n_2}. \qquad (A.7) $$

In our case where the incident wave moves from air to steel, the transmission coefficient is

$$ K_T = \frac{2}{\sqrt{\mu_{r_2}}}, \qquad (A.8) $$

where $\mu_{r_2}$ is the relative permeability of the steel. The transmitted magnetic field strength $H_T$ is therefore

$$ H_T = \frac{2H_c}{\sqrt{\mu_{r_2}}}, \qquad (A.9) $$

where $H_c$ is the magnetic field strength emitted by the coil current. Thus, a value for the $K_1$ constant found in Table 4 can be found using

$$ K_1 = \frac{2}{l_{\text{coil}}\sqrt{\mu_r}}, \qquad (A.10) $$

where it may be used as an initial guess for fitting the parameters to measured data.
### A.0.2 Pseudo Code

In terms of actual code, the following two blocks of simile Python code show how the model is integrated into a simple blocked coil model of the voice coil. 'steelpoleInductor' is the function that calculates the voltage from the influence of the pole. The $H_z$ vector should always be larger than $B_z$ by two values in order to accept the stimulus and a single value due to the boundary conditions of the model. Line 2 provides the stimulus to the system, line 3 converts $B_z$ to $H_z$ using the magnetic permeability and line 4 updates the values of $H_z$ based on the symmetry. Lines 6 and 7 perform two first order spatial gradients, which are then used to estimate the second order spatial gradient on line 9. Line 10 is the $\frac{1}{r}$ first order gradient contribution. Line 12 is the contribution from the point on the axis of symmetry. Line 14 is a simple Euler forward integration scheme. Finally, lines 16 calculates the total time varying magnetic field inside the coil due to the steel pole. The sum in line 16 only extends to N-2 as the value for r is 0 on the axis of symmetry at $n = N - 1$.

0 function u_pole=steelpoleInductor(i1,...
1         K1, K2, K3, mu, Ts, dR, R)
2     Hz(0) = K1*i1
3     Hz(1:N) = Bz./(mu)
4     Hz(N+1) = Hz(N-1)
5
6     dHzDrN = (Hz(1:N) - Hz(0:N-1))
7     dHzDrN1 = (Hz(2:N+1) - Hz(1:N))
8
9     dBzDt(0:N-2) = K2*((dHzDrN1 - dBzDrN)
10                      + dR.*dHzDrN./R(0:N-2))
11
12    dBzDt(N-1) = 4*K2*(Hz(N-1) - Hz(N))
13
14    Bz += Ts*dBzDt
15
16    u_pole = K3*sum(dBzDt(0:N-2).*R(0:N-2))

The next section of code shows the differential equations of the blocked coil model, with a DC resistance $R_e$, stimulated by a voltage *stim*.

0 Fs = 200e3
1 Ts = 1/Fs
2
3 %% System Definitions & Parameters (S.I units)
4 #%Coil Definition
5 nTurns = 78
6 lCoil = 12e-3
7 #Pole Dimensions
8 rPole = 15.9e-3;
9
10 #Spatial discretization
11 nR = 128
12 r = np.linspace(0, rPole, nR)
13 dR = r[1] - r[0]

14 R = (rPole - r)

15

16 #Electrical parameters

17 Re = 2.98

18 Lleak = 3.6e-4

19 LMag = 2.11e-3

20

21 #Material Parameters

22 mu0 = 4*pi*1e-7

23 muRD = 299;

24 sigma = 5.8e6

25

26 #Factors for Method of Lines model

27 K1 = 2/(1Coil*sqrt(muRD))

28 K2 = 1/(sigma*dR^2)

29 K3 = 2*pi*dR

30 mu = mu0*muRD

31 #Sweep definition (User defined here)

32 #stim is the sweep signal of length numSamples.

33 %Solving

34 i = 0

35 i2 = 0

36 for t=0:(numSamples-1)

37

38 i1 = (i - i2)*nTurns

39 u_pole = steelpoleInductor(i1, K1, K2,...

40 K3, mu, Ts, dR, R)

41 didt = (stim(t) - Re*i...

42 - nTurns*u_pole)/Lleak

43

44 di2dt = u_pole*nTurns/LMag

45

46 i += Ts*didt

47 i2 += Ts*di2dt

48 end