A computationally efficient hybrid 2D–3D subwoofer model

A subwoofer generates the lowest frequency range in loudspeaker systems. Subwoofers are used in audio systems for live concerts, movie theatres, home theatres, gaming consoles, cars, etc. During the last decades, numerical simulations have emerged as a cost- and time-efficient complement to traditional experiments in the design process of different products. The aim of this study is to reduce the computational time of simulating the average response for a given subwoofer design. To this end, we propose a hybrid 2D–3D model that reduces the computational time significantly compared to a full 3D model. The hybrid model describes the interaction between different subwoofer components as interacting modules whose acoustic properties can partly be pre-computed. This allows us to efficiently compute the performance of different subwoofer design layouts. The results of the hybrid model are validated against both a lumped element model and a full 3D model over a frequency band of interest. The hybrid model is found to be both accurate and computationally efficient.

The subwoofer is responsible for reproduction of the lowest frequencies in a loudspeaker system. The lower frequency limit for human hearing is around 20 Hz, corresponding to a wavelength of 17 m. This means that the characteristic dimensions of any reasonably sized subwoofer are much smaller than the wavelength. For mathematical modelling, the acoustically small dimensions of a subwoofer justify the use of lumped elements, which are both computationally efficient and conceptually easy to grasp. Lumped elements do have their limitations also for acoustically small components, however. An underlying assumption is that each component (such as the speaker diaphragm, the front chamber, or the ports) of the subwoofer has its own acoustic identity that can be computed separately and that this identity is preserved regardless of the other components. This is not the case if the components are situated so close to each other that their acoustic near-fields interlace. In order to make precise modelling, it is necessary to resort to 3D models of the wave propagation in certain regions. Several commercial software packages that integrate CAD and electroacoustics exist, but while they are convenient for analysing a certain design, computations are time-consuming and therefore impractical if, for example, used in an optimisation loop. Each new layout requires a solution of the governing equation for each frequency of interest. To reduce the computational burden, we propose a hybrid mathematical model that significantly reduces the computational cost of evaluating the subwoofer's performance. The rationale is to pre-compute the properties of fixed parts of the domain so that these parts can be represented by acoustic boundary, or interface, conditions. Then, only part of the domain subject to change is computed in 2D whereas fixed parts of the domain, such as the port to the exterior and the transducer's chamber are computed in 3D. The interface conditions are computed for a surface with a resolution sufficiently high to account also for near field effects.
In this study, we use a fourth order bandpass design 1,2 for the subwoofer. In this design, the transducer is mounted in a sealed back chamber and projects into a ported front chamber, which acts as a passive low pass acoustic filter and only allows a band of frequencies to pass through the output port. The passive acoustic filter has the advantage over an electric filter of being placed so that distortion by the transducer can be filtered out. For more details on bandpass design, readers are referred to Collom's book 3 [p 169] on loudspeaker design. In the following sections, we present three subwoofer models in decreasing order of complexity: 3D, hybrid, and lumped. We validate the results of the hybrid model against the full 3D model and also compare them to a classical lumped element model to assess the hybrid model's efficiency and accuracy.

3D mathematical model
Here, we present a full 3D model that describes the interaction between the subwoofer's different parts. The model uses linear electrical and mechanical circuit models for the transducer, a stiff approximation of the loudspeaker diaphragm, and linear full-wave air acoustics. Such a model is based on fundamental principles of operation for a subwoofer and serves as a stepping stone for all further work. The model thus takes full-wave effects into account but is restricted to the linear working regime, which means that the model will be accurate in the www.nature.com/scientificreports/ small-signal regime. Of nonlinear effects, those related to air acoustics are by no means the most important in a loudspeaker. Instead, displacement-varying stiffness of the diaphragm, magnetic saturation, and the coil leaving the homogeneous region of the magnetic field are examples of nonlinear effects that are more important when leaving the small-signal regime. Another effect not considered here is the modal breakup of the cone. Note that all these complications concern the transducer only and not the 2D/3D effects of air acoustics, which is the main concern here. We consider a subwoofer box standing on an infinite floor as illustrated in Fig. 1. The subwoofer box has dimensions w × h × d , where a transducer is mounted in a sealed back chamber, and an output port of dimensions l × d is placed in the front chamber. The downward facing transducer is mounted at the centre of a baffle with dimensions w × d.
The back chamber, the front chamber, and the exterior of the subwoofer, as illustrated in Fig. 1b, are denoted by b , f , and e , respectively, and b ∪ f is denoted by . Let Ŵ p be the boundary that corresponds to the output port, Ŵ s be all wall boundaries, and Ŵ c be the boundary corresponding to the transducer's cone. Let e c be a unit vector in the axial direction of the speaker pointing into f . The front and back side of the speaker cone are denoted Ŵ c f and Ŵ c b , respectively. We consider time harmonic wave propagation in the linear regime, where the acoustic pressure is given by P(x, t) = R e iωt p(x) , in which ω denotes the angular frequency, t the time, and p the complex pressure amplitude. Under these assumptions, p satisfies Helmholtz equation in ∪ e where c is the speed of sound, k = ω/c is the wave number, and = ∇ · ∇ is the Laplace operator.
We assume that there is no external source and that the complex amplitude p satisfies the Sommerfeld radiation condition, that is uniformly in all directions.
The acoustic pressure is related to the acoustic velocity field u through the linearised Euler equation where ρ is the static air density. Evaluating Eq. (3) at a point on the boundary ∂� with outward directed normal n yields Expression (4) can be used to assign suitable boundary and interface conditions. Assuming all walls to be sound-hard, that is, equipped with boundary condition n · u = 0 , Eq. (4) implies the condition We assume a stiff cone with motion only in the e c direction illustrated in Fig. 1b. At the cone boundary, the normal component of the acoustic velocity will agree with the normal component of the cone velocity. Equation (4) then yields conditions (1) k 2 p + p = 0,   www.nature.com/scientificreports/ where p b and p f are the complex pressure amplitude in b and f , respectively, and n b and n f = −n b are outward directed unit normals on Ŵ c with respect to b and f , respectively. Note that u c is the velocity amplitude, and the cone velocity is u c e c . Next, we consider a lumped mechanical model of the transducer. By Newton's second law, the displacement amplitude δ c of the cone satisfies where M md is the moving mass of the cone, R ms is the cone's mechanical resistance, C ms is the cone's mechanical compliance, and Bl and I denote the force factor and the current, respectively. Multiplying equation (7) by iω and using that u c = iωδ c , yields To close the system, the electromechanical coupling is modelled using the simple electric circuit illustrated in Fig. 1c, assuming that the amplifier delivers a constant voltage V. That is, where R and L denote the resistance and inductance of the voice coil.
Discretisation. Port treatment. Helmholtz equation (1) holds inside the subwoofer as well as in the infinite exterior domain e . However, for computational efficiency, we truncate the domain, compute only in the bounded domain , and represent the response from the exterior through a boundary condition on Ŵ p . For this, we split the boundary Fig. 1b and assume that the acoustic velocity is constant on each panel. Physically this assumption corresponds to a port divided into N p stiff and mass-less pistons. The exterior response is represented through a square impedance matrix Z p 3D of order N p that maps the piece-wise constant acoustic velocities to the mean acoustic pressure on each panel. That is, j is the velocity amplitude in direction n p on Ŵ p j , and p Ŵ p j is the average pressure on Ŵ p j . For the numerical experiments in this paper, we use piecewise constant boundary elements to approximate the problem. More precisely, to assemble the impedance matrix Z p 3D , we solve N p exterior problems by the boundary element method 4 . The algorithm for calculating the impedance matrix, column by column, is as follows: For j = 1, 2, . . . , N p • Set u p to the jth unit vector, this corresponds to setting n · u = 1 on Ŵ p j and n · u = 0 on ∂�\Ŵ j . • Solve Helmholtz equation (1) in e with boundary condition (4) on ∂� and the Sommerfeld radiation condition (2). • Evaluate the average pressure on Ŵ p i , for i = 1, . . . , N p ; the vector of average pressures Finite element discretisation of the subwoofer's interior. Multiplying Helmholtz equation (1) by test functions q f and q b , integrating over f and b , respectively, applying integration by parts, and using the boundary conditions (4), (5), and (6), we obtain and We will apply the finite element method (FEM) to solve discretised versions of the problems (11) and (12) in domains f and b , respectively. Note that problems (11) and (12) are coupled through the cone velocity on Ŵ c .
We triangulate f and b , respectively, and remark that the meshes on f and b do not need to match on their common interface. For the numerical experiments in this paper, we use second order Lagrange basis functions, also known as 10-node tetrahedron elements. We approximate the complex amplitude function of pressure p f and p b , and the test functions q f and q b by www.nature.com/scientificreports/ where N f and N b are the number of degrees of freedom for the finite element approximation in f and b , respectively. In the expressions above θ 1 , θ 2 , . . . , θ N f denote the basis functions in f , and ϕ 1 , ϕ 2 , . . . , ϕ N b denote the basis functions in b .
To describe the subwoofer system, we need to determine , and u p , u c , and I are defined as above. The discretised form of Eqs. (11) and (12) are and Equation (14) is equivalent to To couple p f and u p , we use Eq. (10), which can be written The discretised form of mechanical equation (8) which models the speaker cone displacement can be written where the 1 × N f vector a cf , and the 1 × N b vector a cb have entries To close the system, we need the electric circuit equation (9). Altogether, the full system in matrix form is www.nature.com/scientificreports/

Hybrid mathematical model
In the hybrid model, we split the interior of the subwoofer into two domains (upper and lower) by a horizontal plane at the baffle. Let Ŵ d correspond to the boundary between these two domains. Figure 2a illustrates that the upper box, whose air region is denoted by , is the region above the baffle and contains the transducer, and that the lower box, denoted by , contains the internal air region below the baffle. Note that any internal walls inside the lower box are not included in ; also, we assume that any such walls are fully extruded in the z-direction. In other words, the material layout is planar symmetric along the z-axis. However, in the upper box, the material layout will not be planar-symmetric. To save computational time, we model the wave propagation within the lower box in 2D. The use of the 2D model for the lower box is justified due to planar symmetry, because the frequencies of interest correspond to wavelengths that are long compared to the lower box, and because the waves sufficiently far away from the cone do not exhibit 3D characteristics. However, we rely on 3D models for the upper box and the exterior since 3D effects cannot be avoided there. The upper box, the lower box, and the exterior are treated as separate modules that interact with each other only through interface conditions. Acoustic properties of the 3D models are pre-computed in the form of impedance matrices, Z d and Z p 2D for the upper box and the exterior, respectively. After we have computed the impedance matrices, the lower box is our computational domain where we solve the system of equations using 2D FEM. Computations for various 2D layouts of material in the lower box will then be very fast.  www.nature.com/scientificreports/ and where u d j is the velocity amplitude in the direction of n d on Ŵ dm j .

The mechanical equation (8) for the cone dynamics in its discretised form is
The coupling between p and u on the boundary Ŵ dm depends on the transducer dynamics that is located inside the upper box. Therefore, the impedance matrix also takes into account the cone velocity and current I. The boundary Ŵ dm is divided into N d depth running strips illustrated in Fig. 3a. Similarly, as for the port, we assume constant acoustic velocity on each strip, which physically corresponds to a boundary Ŵ dm made up of stiff massless pistons. The effects of the upper box are then taken into account by an acoustic impedance relation is the average pressure on Ŵ dm j . The matrix Z d maps the N d + 1 vector of panel velocities and cone velocity to the N d + 1 vector of average pressures on the panels and the voice (31) www.nature.com/scientificreports/ coil current. The variables u d i , i = 1, 2, . . . , N d , and u c are the input data for the computations. The impedance matrix Z d is computed in an analogous manner as the algorithm presented for the port treatment on page 3.

Lower box.
We now present a model for , the interior region of the subwoofer's lower box with dimensions w × h l × d illustrated in Fig. 2a. The governing equations are Helmholtz equation (1) in for the acoustic pressure, coupling conditions (4) to the upper box on Ŵ dm and to the exterior on Ŵ p , and sound hard (5) boundary condition on Ŵ s .

The variational form of the governing equation to compute p in is
Here, we assume that acoustic pressure p is planar symmetric in the z-direction. In Fig. 2b, ψ denotes the vertical cross section of , and similarly γ p and γ d denote the vertical cross section of Ŵ p and Ŵ dm in 2D, respectively. The coupling with the upper box is through the impedance matrix Z d , computed by algorithm presented in the previous section, and the coupling with the output port is through the impedance matrix Z Under the above assumptions, the variational form for 3D problem in is reduced to the following 2D problem of finding p ∈ H 1 (ψ) such that We triangulate ψ and define Q basis functions ϑ 1 , ϑ 2 , . . . , ϑ Q . For the numerical experiments in this paper, we use standard continuous, elementwise bi-quadratic basis functions. Here, we determine By putting all the above together, we obtain the following system of governing equations respectively.

Lumped element model
We present a fully lumped element model of the subwoofer modelled by the mass-spring-damper system and the electromechanical circuit, illustrated in Fig. 4. We denote the volume of the back and front chamber by V b and V f , respectively. Let S c and S p be the area, δ c and δ p be the displacement, and u c = iωδ c and u p = iωδ p be the velocities associated with the speaker cone and the output port, respectively. The mechanical force F m , the acoustic force S c p , and the electrical force BlI comprise the total force F t = ma = iωM md u c on the speaker cone membrane: www.nature.com/scientificreports/ Here, a mass spring damper system models the mechanical force is the sum of stiffness due to the speaker and the air in the back chamber. Now, by substituting the current I (9) and the mechanical force F m (38) in Eq. (37), we obtain where the speaker cone impedance is Z c = iωM md + R ms + K t /iω , and the electrical impedance is Z e = R + iωL .
We assume an isentropic flow using the ideal gas law, which states that total pressure times volume raised to power γ (heat capacity ratio) is constant, which entails that where p is the acoustic pressure, and dV f is the change in volume of the front chamber due to movement of speaker cone.
For the lumped model, we assume that u

Validation and discussion
To validate the results of the hybrid model, we consider two layouts of the subwoofer, illustrated in Fig. 5, with width 80 cm, height 100 cm, and depth 60 cm. The downward-facing 18 in. transducer is mounted on a baffle with dimensions 80 cm × 60 cm, located 30 cm from the top. The subwoofer layout on the right in Fig. 5 has a 70 cm wide and 2 cm high horizontal wall that is 48 cm from the floor. We use a generic transducer with realistic parameters, as given in Table 1. The input voltage amplitude to the system is V = 1 V.
The sound pressure level (SPL), 1m in front of the output port, is calculated to assess the performance of the subwoofer. The SPL uses a logarithmic scale to measure the effective pressure of sound relative to a reference pressure in units of decibels (dB). First, we calculate the pressure at 1m by expression (44) p 1m = g T u p , www.nature.com/scientificreports/ where p 1m is the pressure at the floor, 1m in front of the output port, and the N p × 1 vector g is a vector that is pre-computed by the boundary element method while assembling the corresponding port impedance matrix Z p . Now, we can calculate the SPL of our subwoofer design by where p o = 20µPa is the standard reference pressure for SPL calculations. The 3D mathematical model was implemented in COMSOL multi-physics 5 by using unstructured tetrahedral second order elements with maximum edge length h = 0.09 m , the hybrid model was implemented in MATLAB 6 by using uniform square bi-quadratic elements with side length h = 0.01 m , and the lumped model was also implemented in MATLAB. The impedance matrices Z p 3D and Z p 2D for the output port were computed in PAFEC 7 , and the impedance matrix Z d for the driver boundary in COMSOL multi-physics. Figure 6 shows the SPL measured as a function of frequency for the subwoofer layouts in Fig. 5 Voice coil resistance R (�) 5.5 Voice coil inductance L (mH) 1.5 Figure 5. Two subwoofer setups with top-mounted drivers and output ports on the lower right. The subwoofer on the right has an added horizontal wall inside the cabinet. www.nature.com/scientificreports/ closely match the results of the hybrid model for both subwoofer layouts. This validates our assumption of 2D wave propagation within the lower box. For the subwoofer layout, with no material inside the front chamber, the results of the lumped element model are close to the results of the 3D model and follow the same general trend. However, for the second subwoofer layout, with a horizontal wall inside the front chamber, we only compare the results of the 3D and the hybrid model. The lumped model has much more severe limitations compared to the hybrid and the 3D model; lumped elements may work if the elements can be accurately characterised, but doing so may require 2D or 3D calculations in itself for every new subwoofer layout. Therefore, the hybrid and the 3D model are more convenient to use when evaluating a new subwoofer layout. The 3D model accurately evaluates the performance of general internal subwoofer layouts, and the hybrid model accurately and efficiently evaluates the performance of general planar-symmetric subwoofer layouts. Traditionally, the lumped element model has been used to evaluate the performance of subwoofers because it is simple and easy to implement. However, its elements require a new model for every new subwoofer layout, and is thus not suitable to simulate complex layouts. For the 3D model, we only need to pre-compute the impedance matrix Z p 3D to represent the response of the exterior. For the hybrid model, we need to pre-compute the impedance matrices Z p 2D to represent the response of the exterior, but we also need to pre-compute the impedance matrix Z d to represent the response of the upper box. This is a one time cost we pay to compute the impedance matrices in the hybrid model. Then, to evaluate the performance of a given planar-symmetric subwoofer layout, we only need to solve the 2D problem in the lower box. Given that we have calculated the impedance matrices, the hybrid model is significantly faster in computing the wave propagation problem than the full 3D model. For the experiments in this study, we used a system running Windows 10 with Intel Core i7-4790 CPU @ 3.60 GHz processor and 32 GB installed RAM. We performed 20 simulations for time analysis of the 3D and the hybrid model. For each simulation, we used 21 frequencies in the range 20-120 Hz to evaluate the performance of a given subwoofer layout. When using the 3D model, the minimum time for the simulation was 18 min, the maximum time was 20 min 50 s, and the mean time was 19 min 12 s to evaluate the performance at 21 frequencies for a given subwoofer layout compared to the maximum and the minimum time of 42 and 52.5 s, respectively, and the mean time of 44.1 s for the hybrid model. The difference in computational time is significant between the 3D and the hybrid model when we have to evaluate the performance of large sets of subwoofer layouts at multiple frequencies instead of a single frequency. The memory usage for the 3D model is also higher, as it requires 1.25 GB of memory compared to 0.2 GB for the hybrid model. Moreover, we can easily solve finer mesh resolutions for the hybrid model by using elements with side length of h/2, h/4 , or h/8. However, for the 3D model, the system mentioned above runs out of memory when using elements with edge length of h/4, as degrees of freedom reach approximately 1.3 million. With a significant reduction in computational time, the hybrid model is especially beneficial for use in an optimisation loop or when evaluating different planar-symmetric subwoofer layouts to choose the best performing layout.