Finite Element Time-Dependent Analysis
In Part 1 of this course, we discussed the applicability of frequency-domain and time-dependent analyses. Let's recap the circumstances that require a transient model:
- The excitation is a short pulse or includes more than one frequency component.
- The problem becomes nonlinear. The nonlinearity can be due to large signal excitation and higher harmonic generation, a nonlinear magnetization model, or a geometry/topology change because of large deformation.
- The output of interest is the transient response to an excitation.
In this part of the course, we use the same dynamic speaker driver as in Part 2 (which focuses on a frequency-domain study) to discuss the specific considerations for running a transient analysis. The model is solved in the 2D axisymmetric space dimension, and the focuses of this study are to characterize the nonlinear behavior of the driver due to large deformation of the moving parts and to determine the total harmonic distortion (THD) of the acoustic signals produced by the system.
Geometry
The geometry for the transient analysis is shown below. The layer highlighted in blue is created for applying the Perfectly Matched Layer feature in order to avoid unphysical reflections from the artificial boundaries added to form a finite-sized computational domain. The Perfectly Matched Boundary feature used in Part 2 is only available in a frequency-domain pressure acoustics model. For time-domain simulations, use the Perfectly Matched Layer feature, which is applied to a virtual domain attached to an open boundary.
The geometry for a transient analysis. To apply the Perfectly Matched Layer feature to absorb outgoing waves, having a virtual domain attached to the open boundary is necessary.
The geometrical thickness of the perfectly matched layer (PML) in a geometry needs to be set adequately. The thickness can be determined by the following two meshing considerations:
- Use the same mesh element size as that in the adjacent physical domain
- Use a structured mesh with at least 6 mesh layers for the rational scaling and 8 for the polynomial scaling
The recommendation of using the same mesh size in both physical and PML regions is because of the implementation of PML in the time domain: There is no stretching with respect to the real coordinate in the PML, thus, the waves propagate through the PML with the same wavelength that they have in the physical region.
Required Physics Interfaces
Performing the transient analysis requires the Magnetic Fields interface in the AC/DC Module and the Acoustic-Solid Interaction, Transient multiphysics interface in the Acoustics Module. Note that COMSOL® has distinct physics interfaces that solve specific formulations of the wave equation in fluids for different types of analyses. It is therefore important to ensure that the correct physics interface is selected to match the study type. The following image shows the location of the transient Acoustic-Solid Interaction interface in the Select Physics window of the Model Wizard.
The Select Physics window with the Acoustic-Solid Interaction, Transient interface selected and a description of the feature.
The Acoustic-Solid Interaction, Transient multiphysics interface is under the Acoustics > Acoustic-Structure Interaction branch of the Select Physics window of the Model Wizard.
This feature combines the Pressure Acoustics, Transient interface, a transient acoustic model, and the Solid Mechanics interface. Similar to how we performed a frequency-domain analysis in Part 2, the Acoustic-Structure Boundary coupling feature is automatically added to the model to connect the Pressure Acoustics, Transient interface with the Solid Mechanics interface. Later, a Lorentz Coupling feature should be manually added to the Multiphysics node in the model tree to couple the Magnetic Fields interface and Solid Mechanics interface (in the same way as described in Part 2).
Using a Moving Mesh in a Transient Model
The Lorentz Coupling feature is formulated differently for a transient analysis than that for a frequency-domain analysis. In the frequency domain, the formulation explicitly includes the back EMF contributions: It passes the Lorentz force from the Magnetic Fields interface to the Solid Mechanics interface and passes the structural velocity that leads to the induced electric field from the Solid Mechanics interface to the Magnetic Fields interface. However, in the time domain, only the Lorentz force is passed from the Magnetic Fields interface to the Solid Mechanics interface, but the induced currents are not included in the formulation. You need to use a moving mesh in order to capture the back EMF contributions. This specific handling of transient analysis is to avoid double counting of the back EMF when one wants to use a moving mesh to capture nonlinearity due to large deformations and geometry/topology changes, like when the coil is moving partially or fully in and out of the magnetic gap. Since the induced currents are captured automatically by a moving mesh, including them in the coupling formulation will give wrong results.
Therefore, in transient studies it is always necessary to use a moving mesh even for small displacements, because the induced currents are only captured that way. In the current case, we are also interested in modeling the effect of structure deformations on the magnetic field and acoustic field, which adds another source of nonlinearity to the problem.
The Moving Mesh feature can be added by clicking Moving Mesh and choosing Domains > Deforming Domain in the Definitions toolbar, as seen in the image below.
The top half of the COMSOL UI, with the Moving Mesh feature selected in the toolbar and the resulting options with Deforming Domain selected.
The Moving Mesh feature can be added from the Definitions toolbar.
This adds a Moving Mesh interface to the model tree. The following steps are necessary in order to set up the moving mesh correctly:
- In the Deforming Domain 1 subnode, select the air domains above and below the moving parts for Domain Selection and choose Laplace from the Mesh smoothing type list for Smoothing, as seen in the image below. Once the deforming domains are specified, prescribed mesh displacements using the solid displacement components are automatically added to all boundaries next to a solid, and a fixed boundary is used for all other outer boundaries.
The COMSOL Multiphysics UI showing the Model Builder with the Deforming Domain feature selected and the corresponding Settings window with Laplace selected for the mesh smoothing type.
The Laplace mesh smoothing type is used in the deforming domains.
Compared to other nonlinear smoothing methods, the Laplace smoothing is the cheapest option in terms of computations because it is linear and uses one equation for each coordinate direction, which are not coupled together. This method is suitable when the mesh deformation is not too large, and especially when it is used together with an Automatic Remeshing feature, which ensures a new mesh is created whenever the mesh quality becomes less than a given limit.
- Add a Symmetry/Roller feature to specify a symmetry condition on boundaries at r = 0. This allows the mesh to "slide" there.
The COMSOL Multiphysics UI with the Symmetry/Roller feature selected and the corresponding settings and model.
The Symmetry/Roller condition is applied to boundaries lying at r = 0.
Time-Domain Perfectly Matched Layer Settings
The recommended values of the PML scaling factor and scaling curvature parameter are 1,3 for the Polynomial stretching type, and 1, 1 for the Rational stretching type. Here, we use the Polynomial option, as shown in the image below.
The COMSOL Multiphysics UI showing the Model Builder with the Perfectly Matched Layer feature selected, the corresponding settings, and the geometry with the outer layer highlighted in blue to show where the PML feature is applied.
The Perfectly Matched Layer settings for the transient analysis.
Input Signal for Transient Analysis
Performing a transient analysis of the speaker requires a time-dependent voltage signal, V(t), applied to the voice coil. The choice of a specific type of input signal depends on the characteristic under investigation. In this case, the input is a harmonic voltage signal:
The input signal frequency and the amplitude are chosen as f0=70Hz and V0=10V, respectively.
When solving a transient model, it is good practice to ramp up the load to minimize numerical error. An instantaneously applied load leads to a lot of high-frequency content generation due to the inertia effect of the system it applies to, and when they are unresolved, numerical error or even solver divergence can occur. This is also physically motivated since most loads cannot truly change instantaneously in time. You can use a Step or Ramp function with appropriate smoothing for this purpose.
In this analysis, we use the driving signal period, T0=1/f0, to define a ramp function, rm1(t), that has a slope of 1/T0. The ramped input signal, V0*sin(2πf0t)*rm1(t), is then used to define the coil voltage, as shown in the image below. This allows the load to take one period to reach its full magnitude.
The Ramp function selected in the Model Builder, the corresponding Parameters and Smoothing settings, and the function plot.
A close-up of the Coil feature selected in the Model Builder and the Coil section of the corresponding Settings window, with a close-up of the model in the Graphics window.
A Ramp function is defined using the driving signal period (top), and is then used to smooth the input signal (bottom).
Time Stepping Control Setup in Physics Settings
When solving transient models, you should first determine the maximum frequency you want to resolve (known as fmax). Enter the value of the maximum frequency into the Maximum frequency to resolve field in the Transient Solver and Mesh Settings section of the Pressure Acoustics, Transient main physics Settings window, as shown in the image below. The Time stepping setting can also be found in this section, enabling you to set the time-stepping method as either Fixed (preferred) or Free. In this case, it is recommended that you use the Fixed (preferred) method, which is best suited for handling wave propagation. With the settings shown here, the autogenerated solver will be sufficient in most scenarios provided that the computational mesh also resolves the model's frequency content. Here, we set fmax to 3f0 to capture up to the 3rd harmonic.
A close-up of the Pressure Acoustics, Transient feature selected in the Model Builder and the corresponding Settings window showing the Maximum frequency to resolve and Time stepping settings.
The transient solver settings at the top physics level.
Study Steps
Similar to the frequency-domain analysis, the transient study is also carried out using two study steps. First, a Stationary step is used to solve only the electromagnetic part of the problem in order to evaluate the field of the permanent magnet, with the driver in standstill. Then, the full Time Dependent study step takes care of all the relevant multiphysics interactions of the moving speaker.
The Stationary Study selected in the Model Builder and the corresponding Settings window showing the tables in the Physics and Variables Selection setting.
The Time Dependent study selected in the Model Builder and the corresponding Settings window showing settings such as the output times and physics and variables selection.
Study step settings for transient analysis.
The Time Dependent study step solves for the period of time from 0 to T_end = 4T0. The full output time range is split into two subintervals: from 0 to 3T0 and from 3T0 to T_end. The second subinterval is of the most interest since it is assumed that the steady state takes place from time 3T0. Therefore, the output time stepping is finer over the second subinterval in order to obtain a steady-state solution with better time resolution.
The Automatic remeshing option is enabled in the Time Dependent study step to avoid using highly distorted meshes, which could lead to numerically ill-posed problems. In this model, a new mesh is created as soon as the distortion of the current mesh elements exceeds a certain level, which you specify in advance.
A close-up of the Time Dependent study selected in the Model Builder and the Study Extensions settings in the corresponding Settings window.
The Automatic remeshing option is enabled in the Study Extensions section of the Time Dependent study Settings window.
Here are the steps for specifying the condition for remeshing:
- Click Show Default Solver in the Study toolbar.
- Expand the solver sequence and click Automatic Remeshing under the Time-Dependent Solver 1 node.
- In the Condition for Remeshing section of the Settings window for Automatic Remeshing, choose Distortion for the Condition type and enter 2.5 in the Stop when distortion exceeds text field, as seen in the image below.
The Automatic Remeshing feature selected in the Model Builder and the corresponding Settings window showing the Condition for Remeshing settings, such as Distortion set as the condition type.
The Condition for Remeshing setup.
This condition ensures that the solver remeshes when the square root of the maximum element distortion, stored in the variable named comp1.spatial.I1isoMax, becomes larger than 2.5.
As seen in the image below, in the Time Stepping settings, the default transient solver uses the input value of the Maximum frequency to resolve setting to resolve up to the maximum frequency component.
The Time-Dependent Solver selected in the Model Builder and a close-up of the Time Stepping settings in the corresponding Settings window.
Time stepping used by the autogenerated solver.
Initial Mesh
We need an initial mesh to start with, and a new mesh will be automatically generated whenever the mesh distortion tolerance is reached.
For transient models, the mesh needs to be able to resolve the maximal frequency you want to capture, say fmax. This frequency translates to a minimal wavelength, λmin = c/fmax, and in turn to a maximum element size, hmax < λmin/5. Use this criteria to determine the maximum element size for each domain.
The initial mesh is seen in the image below. It uses Mapped meshes for the thin solid domains, including the coil, former, spider, cone, and suspension, and Free Triangular meshes for the rest of the physical domains. Mapped meshes with a distribution of 8 elements are applied to discretize the PML domains. Refer to the attached MPH-file for the detailed meshing sequence.
A close-up of the initial mesh and the mesh of the half-circle.
The initial mesh used for the transient analysis.
The present model runs at a relatively low frequency (including the harmonics) and does not require a resolution of eddy currents (the skin depth) in the pole piece. Thus, compared to the frequency-domain analysis, the boundary layer mesh has been removed.
Transient Solutions and Analysis
The animation below shows the magnetic field in and around the voice coil gap, the motion of the voice coil (orange), the former, the spider, and the cone (pink), as well as how the moving mesh adapts the mesh to the changed topology of the system during the time period of 3T0 to 4T0.
Magnetic field in and around the voice coil gap and the motion of the moving parts and mesh deformation during the time period of 3T0 to 4T0.
The plot below depicts the relative position of the voice coil, calculated as the vertical displacement from its original position averaged over the coil domain. The voice coil is at the lowermost position at 0.044 s (the green marker). At 0.048 s (the red marker) and 0.052 s (the cyan marker), the coil is on its way up and down, respectively.
Relative position of the voice coil.
The distribution of the acoustic pressure around the speaker driver at these three time steps — 0.044 s (left), 0.048 s (middle), and 0.052 s (right) — are shown in the figure below. You can also see how the PML attenuates the acoustic waves generated by the motion of the speaker cone.
Three side-by-side images showing the acoustic pressure at different time steps.
Acoustic pressure at 3 different time steps.
Total Harmonic Distortion (THD) Calculation
An important output of interest of this analysis is the total harmonic distortion of the acoustic signals produced by the system. Distortion analyses, including the total harmonic distortion (THD) and intermodulation distortion (IMD), are commonly used to judge the quality of a speaker driver. These nonlinear effects produce overtones and additional nonharmonic frequency components, the extent of which can be quantified and calculated using the approaches discussed in our blog post "How to Perform a Nonlinear Distortion Analysis of a Loudspeaker Driver". To calculate THD, we first need to choose a point of interest as a listening point, and obtain the frequency spectrum of the steady-state acoustic pressure at this listening point. The value of THD is then calculated as
,
where H{sub}N is the harmonic response of the Nth harmonic and H1 is the fundamental response.
Here, we choose the listening point at r = 0 and z = 0.115 m, the highest point in the physical domain lying on the acoustic axis. It is assumed that the speaker reaches the steady state by time t=3T0, and therefore the pressure signal as a function of time from t=3T0 to t=4T0 at this specific listening point is the steady-state solution that can be used for THD calculation, as seen in the plot below. You can see that the profile slightly differs from a perfect sinusoidal. This means that the output signal also contains frequency components other than f0.
Acoustic pressure at the listening point.
Note that the accuracy of the FFT computation becomes higher for longer time intervals. Therefore, in order to obtain more accurate frequency spectrum of this one cycle signal, we perform a periodic extension of the data over time. This is done by the following steps:
- Use a Point Evaluation feature to evaluate the pressure at the listening point from 3T0 to 4T0 with a time step of T0/50. This creates a table that stores the acoustic pressure as a function of time. See the Pressure at Point node under Results > Derived Values in the attached model file.
- Define an Interpolation function using the Result table as the data source and the table created from Step 1. The signal is interpolated by a cubic spline. See the Interpolation 1 (p_point) node under Global Definitions.
- Define an Analytic function using the interpolation function specified in Step 2 and make it periodic through the Periodic Extension. See the Analytic 1 (p_periodic) node under Global Definitions. The extended pressure signal is seen below.
A periodic extension of the pressure signal at the listening point.
The frequency spectrum analysis of the extended signal is then taken care of by the Global ODEs and DAEs interface defined in Component 2 and the second study, as seen in the images below. The study here consists of a Time Dependent and Time to Frequency FFT step. The first study step solves a global equation that passes the values of the periodically extended signal to a dependent variable named P; the second study computes the FFT of P. The first step is necessary because a Time to Frequency FFT study performs a forward FFT on only a dependent variable that an interface solves for. Here, we choose to extend the signal to 10 acoustic cycles for more accurate FFT computation.
A close-up of Global Equations selected in the Model Builder and the corresponding Settings window with the Global Equations table for variable P.
The periodic extension of the pressure signal is passed onto a dependent variable, P, solved by a Global ODEs and DAEs interface.
The Model Builder with the Time Dependent study selected and the Settings window showing the Study Settings and Physics and Variables Selection table.
The Model Builder with the Time to Frequency FFT study selected and the Settings window showing the Study Settings and the Physics and Variables Selection table.
Study steps used for the frequency spectrum analysis of the extended signal.
The figure below shows the frequency spectrum of the acoustic signal obtained from the FFT study as the value of SPL over frequency. The highest peaks appear at frequencies that are multiples of f0 (red dots), and it is notable that the SPL for the odd-order harmonics is higher than that for the even-order harmonics, which means the symmetrical system nonlinearities dominate the asymmetrical ones.
A graph showing the 1/3 octave bands (PSD) as green bar graphs and the continuous PSD as a blue line, with both lines peaking after 50 Hz.
The frequency spectrum of the acoustic signal (represented as SPL over frequency) at the listening point.
The THD value is computed in Derived Values > THD Evaluation using Global Evaluation with the following expression:
sqrt(sum(with(11 + 10*k, abs(comp2.P)^2), k, 1, 9))/with(11, abs(comp2.P))
The Model Builder with Global Evaluation selected and the corresponding Settings window showing the dataset data and the calculated expression.
THD calculation using Global Evaluation.
The expression uses the with operator to access the solution of a specific frequency component and the sum operator to compute the sum of the expressions for all index values. More specifically, with(11, abs(comp2.P)) provides the acoustic pressure amplitude for the fundamental frequency (the 11th output frequency is 70 Hz), and sum(with(11 + 10*k, abs(comp2.P)^2), k, 1, 9), where k is the summation index, calculates the sum of the squared pressure amplitude of the first 9 harmonics. If you take a look at the output frequencies from the FFT study, you will see the index corresponding to the fundamental frequency and its harmonics. Such a calculation of including 9 harmonics yields a THD equal to 2.5%.
Further Learning
This tutorial model and detailed documentation with step-by-step analysis instructions are also available for download in the Application Gallery: Loudspeaker Driver — Transient Analysis.
Submit feedback about this page or contact support here.