PDF

Helmholtz Resonator Analyzed with Different Frequency Domain Solvers
Introduction
This example simulates a simple three-dimensional axisymmetric Helmholtz resonator, a classic acoustics model of a resonating circuit with a known theoretical solution. The idealized version considered here consists of a tube and a closed volume in series, which are exposed to a harmonically oscillating pressure. Real-world phenomena explained by the resonator include, among others, the resonance from blowing across the top of an empty bottle and the sound produced by closed-cavity drums such as the djembe and by subwoofers. This model illustrates the use of different numerical solvers. The model shows how to solve this pressure acoustics problem for a range of frequencies using the following solvers: (a) Frequency Domain, with and without Asymptotic Waveform Evaluation (AWE) for faster solution, and (b) Frequency Domain-Modal, which reconstructs the frequency response based on a specified set of eigenmodes.
Figure 1: Illustrations of the Helmholtz resonator.
Model Definition
The model consists of a tube and volume coupled in series driven by a harmonically oscillating pressure p0( t) = p0 eiωt at the tube inlet. The Helmholtz resonator (schematically depicted in Figure 1) is one of the simplest resonating circuits. This circuit is typically described using lumped-parameter (equivalent circuit) modeling (see, for example, Ref. 1) as a serial coupling of an acoustic inertance L (equivalent to inductance in electric circuits and to mass in point-mass mechanics) caused by acceleration of the fluid in the tube and an acoustic compliance C (equivalent to capacitance in electric circuits and a spring in point-mass mechanics) arising from compression of the volume; see Figure 1 and Figure 3.
Figure 2: Diagram of the equivalent electric circuit to the Helmholtz resonator.
Because compression of the fluid volume and acceleration of the fluid in the tube are not instantaneous events but instead occur on specific time scales (each given by the geometry together with the properties of the fluid), the response of the resonator depends on the frequency, and there exists a frequency that maximizes the response — in other words, a resonance frequency ωR. This is of course to be expected given the direct analogy to electric LC circuits, and by the same analogy we find that the resonance frequency is given by
(1)
As detailed in Ref. 1, the acoustic lumped-parameter elements are given by
(2)
where ρ0 is the background quiescent density of the fluid, c0 is the background quiescent speed of sound, l is the tube length, a is the tube radius, S is the cross-sectional area of the tube transverse to the direction of the flow, γ is the end correction factor (detailed below), and V is the closed volume. Thus, using Equation 2 in Equation 1 we find that the resonance frequency ωR is given by
(3)
where in the last equation we have assumed that the volume is a sphere with and the tube is cylindrical so S = πa2.
End correction factor γ: physical origin and approximate numerical value
When the fluid exits the tube and enters the volume, the acoustic waves disperse and the acoustic pressure drops. However, the waves initially continue along the axis of the tube when they just leave it, and moreover, they cannot move into the region occupied by the tube. Consequently, they do not completely disperse immediately as they leave the tube and the immediate region downstream of the tube is therefore still felt by the fluid in the tube where it imposes an acoustic load. In ideal models, this load results in an additional acoustic inertance corresponding to an effective increase in the length of the tube by γa. In other words, the total length of the lumped-parameter inertance L in Equation 2 is longer than the actual tube length l.
The factor γ depends on the specific geometry of the tube-volume connection and is of the order unity, Ref. 1 and 2. We shall take for reference an infinite flange for which the correction factor is γ = 0.82 (Ref. 2); this is not fully correct because the acoustic pressure disperses less in the circular geometry resulting in a larger acoustic load.
Solver descriptions
Using the default settings for the Frequency Domain solver, it solves the problem subjected to harmonic excitation at a set of specified excitation frequencies. While this can be time-consuming for larger frequency sweeps, the (numerically) exact solution is calculated explicitly at every frequency, and so the solutions from this solver can always be expected to be correct (assuming convergence of the model and appropriate meshing to resolve all length scales of the physics).
To accommodate larger sweeps, this solver also contains the option to apply asymptotic waveform evaluation (AWE). This approach does not explicitly compute the exact solution at all frequencies but instead performs a Taylor expansion of the solution about a few exact solutions and otherwise uses a lower-order approximation (Padé or Taylor) to estimate the solution across the required frequency range.
Finally, the Frequency Domain Modal solver can also be used to perform a frequency sweep. For this, it first computes a set of system eigenfrequencies and associated eigenmodes (searching either within a user-defined range, or for a user-defined number of frequencies). The full solution across the frequency sweep is then approximated by a linear combination of the basis set formed by the eigensolutions (see Notes About the COMSOL Implementation).
Goals of analysis
This pressure acoustics problem is solved in a specified frequency regime using different solvers, with the dual purpose of illustrating the capabilities of the solvers and also highlight the solver-specific settings to be aware of. The following preset solvers are used: (a) Frequency Domain, with and without Asymptotic Waveform Evaluation (AWE) for faster solution, and (b) Frequency Domain-Modal, which reconstructs the frequency response based on the eigenfrequencies in the specified range.
Resistance due to acoustic radiation
As a final note, it should be mentioned that the equivalent circuit diagram in Figure 2 is incorrect. A full equivalent circuit description of the resonator should also include the acoustic radiation resistance R caused by the dissipated energy out of the tube when the fluid in the tube moves into the volume, see Figure 3. Mathematical formulas for the resistance in a number of situations may be found in Ref. 1. However, the value of the resistance does not affect the resonance frequency (the only read-out used here), only the absolute level of the impedance of the system, so we are well-justified in ignoring R.
Figure 3: Diagram of the equivalent electric circuit to the Helmholtz resonator including the radiation resistance R.
Results and Discussion
This model uses the absolute acoustic pressure averaged over the end volume
(4)
to investigate the response. We find that the lowest eigenfrequency indeed corresponds to the resonance predicted by the lumped-parameter model, and that this theoretical prediction indeed is in good agreement with the numerical results. However, the full numerical COMSOL Multiphysics solutions illustrate the many higher modes that are ignored in the simple lumped-parameter model.
Comparing the different solvers, it is noted that they all produce the same response for the frequency sweep; see Figure 4 and Figure 5 on logarithmic and linear scales, respectively. The sound pressure level at 1043 Hz is depicted in Figure 6. The AWE and modal based frequency sweep methods are quite useful to speed up the running times of large models with a large number of frequencies requested. Due to the small number of degrees of freedom in this model, it is not possible to measure any speed-up in CPU time. See the next section for details.
Figure 4: Frequency sweep for the readout illustrating agreement between all solvers and the theory for the first eigenfrequency and also illustrating agreement between all solvers for the higher frequencies.
Figure 5: Same frequency sweep as in Figure 4 here shown on a linear frequency scale to emphasize the higher frequencies.
Figure 6: Sound pressure level computed as a default output from the simulations, here shown for f = 1043.2 Hz.
Notes About the COMSOL Implementation
When several resonance frequencies are present, the default parameters for the AWE option may use some tuning, for example, using an adequate Absolute tolerance value (see under the >AWE Solver 1 node). In this model, the default value of 0.001 is used and gives reasonable results (when compared to the other methods). Lowering the value to 0.0001 will improve the high frequency results slightly (see the linear scale plot). Other settings are in general less important, like changing the number of points to linearize about (Evaluation points), changing the number of terms in the Taylor expansion about each point (Expansion size), or changing from Padé to Taylor expansion of the approximating solution (Expansion type). In the current setup we chose a relatively high upper frequency bound so about 10 resonance frequencies are present in the sweep. Had the upper frequency limit been lowered so only the first resonance was included (for instance, by setting fmax = 100 Hz), then the default relative AWE tolerance would have sufficed. Note that if either end of the frequency range is close to a resonance, the AWE solver can become unreliable.
The linper operator informs the solver that the term in the expression is a perturbation (a source term) that must be included in the linearized problem. The modal solver will only use the pressure under the linper operator as a source, while the other solvers will ignore this perturbation term.
To obtain good results with the Modal solver up to fmax, we must set the upper limit in its eigenfrequency search to 1.5*fmax to capture modes that may have an influence on the highest part of the frequency of interest. Using only fmax as the upper limit results in poorer estimates of the solution at higher frequencies.
The AWE option and the Frenquency Domain-Modal solver both rely on approximating the solution using a few exact solutions in the sweep range. Thus, these methods provide greater speed-up in CPU time if only few resonances fall in the sweep range, or if it is comparatively easier to find the eigensolutions relative to all the full solutions. Thus, these methods will be especially useful for large models and fine frequency sweeps over broad ranges. As an example, it is possible to update the parameter hmax to 0.5 mm to increase the size of the model. Solving the frequency sweep will take then approximately 170 s, while the AWE solver will require 55 s and the modal based frequency sweep will require only 39 s. The speed up will depend on the number of frequencies requested and the number of resonances present in the frequency of interest.
Finally, note that we define the model in terms of geometric parameters (a, L, and R). This makes it easy to quickly include parametric sweeps in the geometry, which, for instance, could be used to tune the lowest eigenfrequency.
References
1. D.T. Blackstock, Fundamentals of Physical Acoustics, John Wiley & Sons, 2000.
2. A.D. Pierce, Acoustics: An introduction to its physical principles and applications, Acoustical Society of America, 1989.
Application Library path: Acoustics_Module/Tutorials,_Pressure_Acoustics/helmholtz_resonator_solvers
Modeling Instructions
From the File menu, choose New.
New
In the New window, click  Model Wizard.
Model Wizard
1
In the Model Wizard window, click  2D Axisymmetric.
2
In the Select Physics tree, select Acoustics>Pressure Acoustics>Pressure Acoustics, Frequency Domain (acpr).
3
Click Add.
4
Click  Study.
5
In the Select Study tree, select General Studies>Frequency Domain.
6
Global Definitions
Parameters 1
Load all model parameters from a file; these include geometrical parameters and physical properties of the air.
1
In the Model Builder window, under Global Definitions click Parameters 1.
2
In the Settings window for Parameters, locate the Parameters section.
3
Click  Load from File.
4
Geometry 1
Construct the simple 2D axisymmetric geometry of the resonator from Figure 1 using rectangles and circles.
Rectangle 1 (r1)
1
In the Geometry toolbar, click  Rectangle.
You start by making the tube. It has length L and radius a.
2
In the Settings window for Rectangle, locate the Size and Shape section.
3
In the Width text field, type a.
4
In the Height text field, type L.
5
Locate the Position section. In the z text field, type -L.
Circle 1 (c1)
1
In the Geometry toolbar, click  Circle.
Proceed to make the volume which has radius Rv and is placed after the tube. You ensure an overlap of a/2 knowing that you will merge the tube and volume into one object (the resonator). Note that objects or parts of objects extending into the left-hand side of the symmetry line (r = 0) will be removed, so you will not have to manually remove the left half of the circle.
2
In the Settings window for Circle, locate the Size and Shape section.
3
In the Radius text field, type Rv.
4
Locate the Position section. In the z text field, type -(L+Rv-a/2).
Union 1 (uni1)
1
In the Geometry toolbar, click  Booleans and Partitions and choose Union.
Clean up the remaining parts of the geometry.
2
Click in the Graphics window and then press Ctrl+A to select both objects.
Delete Entities 1 (del1)
1
In the Model Builder window, right-click Geometry 1 and choose Delete Entities.
2
On the object uni1, select Boundaries 1, 2, and 5 only.
3
In the Settings window for Delete Entities, click  Build All Objects.
4
Click the  Zoom Extents button in the Graphics toolbar.
5
In the Geometry toolbar, click  Build All.
Materials
Now add the material. For this simple model it suffices to use tabulated values for the physical parameters of air, which can be defined manually by adding a blank material and assign to it the values rho_air and c_air. Alternatively load air from the Material Library.
My air
1
In the Model Builder window, under Component 1 (comp1) right-click Materials and choose Blank Material.
2
In the Settings window for Material, type My air in the Label text field.
3
Locate the Material Contents section. In the table, enter the following settings:
Next define integrated variables and an integration operation. First define an integration coupling that integrates over the volume V which will be used for evaluating the system response intop1. Then define the average absolute pressure in the volume which will be used as a readout of the system response.
Definitions
Integration 1 (intop1)
1
In the Definitions toolbar, click  Nonlocal Couplings and choose Integration.
2
Variables 1
1
In the Model Builder window, right-click Definitions and choose Variables.
2
In the Settings window for Variables, locate the Variables section.
3
Click  Load from File.
4
Pressure Acoustics, Frequency Domain (acpr)
Pressure Acoustics 1
Next, you define that you will be solving Pressure Acoustics in the entire domain, and you apply the boundary conditions.
1
Click the  Zoom Extents button in the Graphics toolbar.
Pressure 1
1
In the Model Builder window, right-click Pressure Acoustics, Frequency Domain (acpr) and choose Pressure.
2
The applied pressure is a combination of an constant part and a linper operator. The constant part will be used in the Frequency Domain steps and inactive in the rest, while the linper part will be active when running the Frequency Domain, Modal Step and inactive in the rest, as described in Notes About the COMSOL Implementation.
3
In the Settings window for Pressure, locate the Pressure section.
4
In the p0 text field, type 1+linper(1).
Mesh 1
Next, define the mesh to have a maximum size corresponding to one fifth of the smallest wavelength in the system (see Notes About the COMSOL Implementation).
Free Triangular 1
In the Mesh toolbar, click  Free Triangular.
Size
1
In the Model Builder window, click Size.
2
In the Settings window for Size, locate the Element Size section.
3
Click the Custom button.
4
Locate the Element Size Parameters section. In the Maximum element size text field, type hmax.
5
Click  Build All.
Study 1 - Frequency sweep
Next, set up the first of the three solvers (Frequency Domain solver). You already included this solver (this study) when you started this analysis, so you can immediately set it up.
1
In the Model Builder window, click Study 1.
2
In the Settings window for Study, type Study 1 - Frequency sweep in the Label text field.
Step 1: Frequency Domain
1
In the Model Builder window, under Study 1 - Frequency sweep click Step 1: Frequency Domain.
2
In the Settings window for Frequency Domain, locate the Study Settings section.
3
In the Frequencies text field, type 10^{range(log10(fmin),(log10(fmax)-(log10(fmin)))/499,log10(fmax))}.
This command selects 500 frequencies in the range fmin - fmax which will be evenly spaced when shown on a logarithmic axis.
4
In the Home toolbar, click  Compute.
Root
Now use the same solver as before, but this time with the Asymptotic Waveform Evaluation (AWE) option. You perform this as a separate study to be able to compare the two solutions, and therefore you first add a new Frequency Domain study.
Add Study
1
In the Home toolbar, click  Add Study to open the Add Study window.
2
Go to the Add Study window.
3
Find the Studies subsection. In the Select Study tree, select General Studies>Frequency Domain.
4
Click Add Study in the window toolbar.
5
In the Home toolbar, click  Add Study to close the Add Study window.
Study 2
Step 1: Frequency Domain
1
In the Settings window for Frequency Domain, locate the Study Settings section.
2
In the Frequencies text field, type 10^{range(log10(fmin),(log10(fmax)-(log10(fmin)))/499,log10(fmax))}.
3
Click to expand the Study Extensions section. Select the Use asymptotic waveform evaluation check box.
4
This function is used during the AWE algorithm to evaluate its performance. Other functions, including the value in a point somewhere in the system, could be used as well.
5
In the Model Builder window, click Study 2.
6
In the Settings window for Study, locate the Study Settings section.
7
Clear the Generate default plots check box.
8
In the Label text field, type Study 2 - Frequency sweep with AWE.
9
In the Home toolbar, click  Compute.
Add Study
1
In the Home toolbar, click  Add Study to open the Add Study window.
2
Go to the Add Study window.
3
Find the Studies subsection. In the Select Study tree, select Preset Studies for Selected Physics Interfaces>Frequency Domain, Modal.
4
Click Add Study in the window toolbar.
5
In the Home toolbar, click  Add Study to close the Add Study window.
You then set up the third solver; the Frequency Domain Modal study.
Study 3 - Modal solver frequency sweep
1
In the Model Builder window, click Study 3.
2
In the Settings window for Study, locate the Study Settings section.
3
Clear the Generate default plots check box.
4
In the Label text field, type Study 3 - Modal solver frequency sweep.
Step 1: Eigenfrequency
1
In the Model Builder window, under Study 3 - Modal solver frequency sweep click Step 1: Eigenfrequency.
2
In the Settings window for Eigenfrequency, locate the Study Settings section.
3
From the Eigenfrequency search method list, choose Region.
Next, change the upper limit of the real part of the frequencies you want to investigate to make sure the entire frequency range fmin - fmax is appropriately resolved.
4
Find the Search region subsection. In the Largest real part text field, type 1.5*fmax.
Step 2: Frequency Domain, Modal
1
In the Model Builder window, click Step 2: Frequency Domain, Modal.
2
In the Settings window for Frequency Domain, Modal, locate the Study Settings section.
3
In the Frequencies text field, type 10^{range(log10(fmin),(log10(fmax)-(log10(fmin)))/499,log10(fmax))}.
4
In the Home toolbar, click  Compute.
Results
Frequency sweeps - logarithmic scale
1
In the Home toolbar, click  Add Plot Group and choose 1D Plot Group.
Next, you add a plot of the system response across all investigated frequencies.
2
In the Settings window for 1D Plot Group, type Frequency sweeps - logarithmic scale in the Label text field.
Frequency sweep
1
Right-click Frequency sweeps - logarithmic scale and choose Global.
2
In the Settings window for Global, locate the Data section.
3
From the Dataset list, choose Study 1 - Frequency sweep/Solution 1 (sol1).
4
Locate the y-Axis Data section. In the table, enter the following settings:
5
In the Label text field, type Frequency sweep.
Frequency sweep with AWE
1
In the Model Builder window, right-click Frequency sweeps - logarithmic scale and choose Global.
2
In the Settings window for Global, type Frequency sweep with AWE in the Label text field.
3
Locate the Data section. From the Dataset list, choose Study 2 - Frequency sweep with AWE/Solution 2 (sol2).
4
Locate the y-Axis Data section. In the table, enter the following settings:
Modal solver frequency sweep
1
Right-click Frequency sweeps - logarithmic scale and choose Global.
2
In the Settings window for Global, type Modal solver frequency sweep in the Label text field.
3
Locate the y-Axis Data section. In the table, enter the following settings:
4
Locate the Data section. From the Dataset list, choose Study 3 - Modal solver frequency sweep/Solution 3 (sol3).
Theoretical reference
1
Right-click Frequency sweeps - logarithmic scale and choose Global.
2
In the Settings window for Global, type Theoretical reference in the Label text field.
3
Locate the Data section. From the Dataset list, choose Study 1 - Frequency sweep/Solution 1 (sol1).
4
Locate the y-Axis Data section. In the table, enter the following settings:
5
Locate the x-Axis Data section. From the Parameter list, choose Expression.
6
In the Expression text field, type f_theo.
7
Click the  x-Axis Log Scale button in the Graphics toolbar.
8
Click the  y-Axis Log Scale button in the Graphics toolbar.
9
In the Frequency sweeps - logarithmic scale toolbar, click  Plot.
Frequency sweep
1
In the Model Builder window, click Frequency sweep.
2
In the Settings window for Global, click to expand the Coloring and Style section.
Now change plot styles and update legend texts to produce presentation-worthy plots of the system response computed from all the solvers. These will look like Figure 3 and Figure 4.
Frequency sweep with AWE
1
In the Model Builder window, click Frequency sweep with AWE.
2
In the Settings window for Global, locate the Coloring and Style section.
3
Find the Line style subsection. From the Line list, choose None.
4
Find the Line markers subsection. From the Marker list, choose Plus sign.
5
In the Number text field, type 70.
Modal solver frequency sweep
1
In the Model Builder window, click Modal solver frequency sweep.
2
In the Settings window for Global, locate the Coloring and Style section.
3
Find the Line style subsection. From the Line list, choose None.
4
Find the Line markers subsection. From the Marker list, choose Point.
5
In the Number text field, type 72.
Frequency sweeps - logarithmic scale
1
In the Model Builder window, click Frequency sweeps - logarithmic scale.
2
In the Settings window for 1D Plot Group, click to expand the Title section.
3
From the Title type list, choose Label.
4
Locate the Legend section. From the Position list, choose Upper left.
5
Locate the Plot Settings section. Select the x-axis label check box.
6
7
Select the y-axis label check box.
8
9
Click to expand the Title section. In the Frequency sweeps - logarithmic scale toolbar, click  Plot.
Frequency sweep
1
In the Model Builder window, click Frequency sweep.
2
In the Settings window for Global, click to expand the Legends section.
3
From the Legends list, choose Manual.
4
Frequency sweep with AWE
1
In the Model Builder window, click Frequency sweep with AWE.
2
In the Settings window for Global, locate the Legends section.
3
From the Legends list, choose Manual.
4
Modal solver frequency sweep
1
In the Model Builder window, click Modal solver frequency sweep.
2
In the Settings window for Global, locate the Legends section.
3
From the Legends list, choose Manual.
4
Theoretical reference
1
In the Model Builder window, click Theoretical reference.
2
In the Settings window for Global, locate the Legends section.
3
From the Legends list, choose Manual.
4
5
Locate the Coloring and Style section. From the Color list, choose Black.
Frequency sweeps - linear scale
1
In the Model Builder window, right-click Frequency sweeps - logarithmic scale and choose Duplicate.
2
Click the  x-Axis Log Scale button in the Graphics toolbar.
3
In the Settings window for 1D Plot Group, type Frequency sweeps - linear scale in the Label text field.
4
Locate the Legend section. From the Position list, choose Upper right.