Acoustics in the Frequency Domain

Introduction
Acoustics is the field of physics that models sound waves by changes in pressure. Two approaches to model acoustic systems are common: One approach is to model acoustics in the time domain and the other is to model in the frequency domain. This tutorial focuses on the modeling of sound in the frequency domain and makes use of the Helmholtz partial differential equation (PDE) as the model. The acoustic modeling in the frequency domain introduced here will build on concepts introduced in the tutorial Acoustics in the Time Domain which should be read as a first exposure to acoustics modeling. Since the two approaches are related the relations will be presented throughout this tutorial.
The main model of acoustics in the frequency domain is the Helmholtz equation. The Helmholtz PDE is a time independent equation. Because the Helmholtz PDE is a time independent PDE it can be solved more efficiently compared to the time dependent wave equation used for modeling acoustics in the time domain. The Helmholtz equation is, however, only applicable when modeling acoustic systems which have a harmonic time dependency. In other words, an inharmonic sound signals have to be modeled in the time domain, and the benefits of using the Helmholtz equation can not be exploited.
Two types of analysis in the frequency domain are introduced in this tutorial: Time Harmonic Analysis and Eigenfrequency Analysis. Both the time harmonic and the eigenfrequency analysis are based on the Helmholtz PDE model in conjunction with various types of boundary conditions, which are also introduced in this tutorial. The purpose of a time harmonic analysis is to compute the frequency response of an acoustic system over a range of frequencies. An eigenfrequency analysis on the other hand is applied to solve for the eigenmodes and eigenfrequencies of an acoustic system. The actual analysis of a time harmonic model is done with ParametricNDSolve and for an eigenfrequency analysis NDEigensystem is made use of.
The symbols and corresponding units used throughout this tutorial are summarized in the Nomenclature section.
Load the finite element package.
In[22]:=
Click for copyable input
Time Harmonic Analysis
Time harmonic analysis is a branch of acoustics concerned with the frequency response of an acoustic system. A sound signal is referred to as time-harmonic if it can be expressed as a sine function with a specific frequency. For a time harmonic analysis an acoustic system is exposed to several harmonic sound signals over a range of frequencies, and the performance of the device at frequency of interest is analyzed. This type of analysis is important when building frequency-dependent acoustic systems, for example a low-pass filter that is designed to attenuate sound at higher frequencies.
In response to harmonic stimulus the resulting sound pressure field can be shown to be time-harmonic as well [1]. A sound pressure field is said to be time-harmonic if the pressure variation at any spatial position has a sinusoidal time dependence with an angular frequency .
A general expression of a harmonic sound pressure field is written as:
Here denotes the amplitude at the given position, and is an initial phase shift at . In cases where two sound signals have no phase difference are said to be "in phase". If the phase difference happens to be , the two sound signals are said to be "in antiphase".
For analytical convenience, the time-harmonic relation (2) is often expressed in complex form known as the complex exponential representation (CER).
The CER expression (3) can be interpreted as a rotating vector in the complex plane. The following figure illustrates this behavior.

"Complex Exponential Representation"

The rotating vector is known as the complex amplitude function. The amplitude function rotates counterclockwise at a speed of the angular frequency . At any given time the projection of on the real axis represents the instantaneous sound pressure, and the vector length corresponds to the local amplitude.
If we express the amplitude function and its complex conjugate as and , then the local amplitude can be calculated by
When the time harmonic relation (4) is inserted into the wave equation the equation simplifies to a time independent Helmholtz equation. The derivation of the Helmholtz equation from a wave equation will be presented in a later section entitled Derivation of the frequency acoustics model from the time domain model. For now it is important to understand that an unknown sound field can be solved for in the frequency domain by using the angular frequency in the Helmholtz PDE model (5).
The terms and are monopole and dipole sources, respectively.
The computed solution , however, can be easily transformed back into time domain using the time-harmonic relation (6).
Eigenfrequency Analysis
If the monopole source and the dipole source are removed from the Helmholtz PDE model (7), the equation simplifies to the source-free Helmholtz equation:
Equation (8) can be treated as an eigenvalue problem such that , and can be solved with NDEigensystem. Here the differential operator corresponds to the left hand side of (9), and represents the eigenvalue of the eigensystem.
The set of eigenvalues that fulfills the source-free Helmholtz equation gives the corresponding eigenfrequencies by
or
The amplitude function that is paired with each eigenvalue is called eigenmode.
The eigenfrequency is also known as the natural frequency, which determines the resonance of an acoustic system. To illustrate the acoustic resonance we consider an open-ended tube and a closed-ended tube. Both tubes are filled with air, and the length .
By convention, the closed-ended tube denotes a tube with one closed end only.

Since the medium pressure must be equal to the ambient reference pressure at the open end, the sound pressure is fixed at zero such that . At the closed end, however, the sound pressure accumulates and reaches its maximum value since no forward motion is possible. Due to these boundary conditions the tube can only sustain sound waves at certain frequencies, that is, the eigenfrequencies.
At each eigenfrequency a standing wave forms within the tube. The lowest eigenfrequency, which corresponds to eigenmode 1, is called the fundamental frequency. As shown in the acoustics time domain tutorial, compared to a traveling wave with the same amplitude a standing wave needs a weaker sound source. In other words, a sound source excites an acoustic system the most at each eigenfrequency.
Eigenfrequency analysis is therefore an important consideration when designing acoustic systems that utilize (or prevent) resonance, such as musical instruments, acoustic filters and concert halls.
Helmholtz Equation
Introduction to Helmholtz Equation
The behaviour of an acoustics system in the frequency domain is investigated by repetitively solving the Helmholtz PDE for a specific frequency out of a frequency range of interest. The Helmholtz equation (10) is used for modeling a harmonic sound pressure field at a specific angular frequency
The dependent variable in the Helmholtz equation is the sound pressure . The sound pressure wave is propagating in a medium with density at the speed of sound . The sound pressure field is modeled in response to a harmonic sound stimulus at a frequency , which is related to the angular frequency by .
Sound pressure can be understood as the local pressure deviation from an ambient reference pressure: , where denotes the position vector. Terms and represent monopole and dipole sources, respectively. The Source Types section describes these sound sources.
Define the acoustic model.
In[26]:=
Click for copyable input
Various sections in the documentation explain the use of inactive PDE operators. Please refer to Numerical Solution of Partial Differential Equations.
Set up a 1D time independent acoustic model in the frequency domain.
In[29]:=
Click for copyable input
Out[29]=
As shown in the tutorial Acoustics in the Time Domain, the transient acoustic model can be set up in a similar way:
Set up a 1D transient acoustic model in the time domain.
In[30]:=
Click for copyable input
Out[30]=
Note that for the frequency domain acoustic model, the time derivative term has been converted to the frequency dependent term by using the time-harmonic relation (11). The derivation can be found in the following section: Derivation of the frequency acoustics model from the time domain model.
To make use of specific material parameters in the equation we extract relevant data from the ThermodynamicData.
In[31]:=
Click for copyable input
The actual material parameter insertion can then be done in the following way:
In[33]:=
Click for copyable input
Out[33]=
The following 1D example shows a frequency domain acoustic model simulation. In the first step a time harmonic analysis will be performed and in the subsequent step an eigensystem analysis is done of the same acoustic mode. The relation between the two analysis types will then be apparent.
A radiation boundary condition is set on the left end to serve as the sound signal input. By default a sound hard boundary condition is implicitly used at the right end.
In[34]:=
Click for copyable input
Out[34]=
Insert the material parameters into the model.
In[35]:=
Click for copyable input
Solve the PDE repetitively from to with an increment of .
In[36]:=
Click for copyable input
Compute the frequency response at the sample frequencies.
In[39]:=
Click for copyable input
Out[40]=
In a next step an eigenvalue analysis is performed.
Solve for the five smallest eigenvalues with NDEigenvalues.
In[41]:=
Click for copyable input
Calculate the corresponding eigenfrequencies with the relation .
In[42]:=
Click for copyable input
Out[42]=
Inspect the eigenfrequencies and the frequency response spectrum.
In[43]:=
Click for copyable input
Out[43]=
Note that at each eigenfrequency the amplitude response reaches its local maximum.
Derivation of the frequency acoustics model from the time domain model
The Helmholtz equation is derived from the wave equation (12) with harmonic time-dependence. The general wave equation is given as:
Here the terms and are monopole and dipole sources, respectively.
Recall that if a sound pressure field is assumed to be time-harmonic, the pressure variation in time for a particular frequency can be expressed in the complex plane by an amplitude function .
Likewise, the monopole and dipole sources can be expressed by amplitude functions and , respectively.
Taking the second order time derivative of (13) gives
Taking the gradient of (14) yields
Insert (15), (16), (17) into (18), the wave equation becomes
Factor out the common term then the equation simplifies to the time-independent, inhomogeneous Helmholtz equation (19).
Model Parameter Setup
In acoustics simulations the wavelength of a sound wave needs to be resolved by a sufficiently fine mesh in order to get an accurate numerical solution of the governing partial differential equation. Here we create a function to set the MaxCellMeasure for time harmonic waves.
For a time harmonic wave the wavelength is . The default resolution of the mesh is set to 20 elements per , which means that the max cell measure would be twenty times smaller than the sound wavelength . Typically this is sufficient to resolve waves accurately [20]. However, it is always possible to assign other resolution values to meet different accuracy requirements.
Define a function to set the max cell measure.
In[44]:=
Click for copyable input
Define a function named AcousticOptions that gives NDSolve options appropriate for solving acoustic problems.
In[45]:=
Click for copyable input
The following model parameters are used for the examples in this tutorial. These parameters define the simulation domain .
Set up model parameters.
In[47]:=
Click for copyable input
In many examples we will be using a radiation boundary condition to produce a harmonic sound wave, and an absorbing boundary condition to avoid wave reflection.
Define
Click for copyable input
and
Click for copyable input
to be a radiation boundary condition and an absorbing boundary condition, respectively. The subscripts and denote the sides where a boundary condition is enforced.
In[50]:=
Click for copyable input
Source Types
The Helmholtz PDE model (21) contains two types of time harmonic sources: monopole and dipole source . The following sections will demonstrate how these sound sources are set up for modeling in the frequency domain. The physical meaning of sound sources is explained in detail in the acoustics time domain tutorial.

Monopole Sources

The monopole source from the Helmholtz equation (22) models a point source that radiates sound isotropically. An example of an acoustic monopole would be a small sphere whose radius alternately expands and contracts [23].
To make use of a monopole source we have to specify the monopole source strength and the source location . The monopole source term may be written as
where is a regularized Dirac delta function at the source location.
There are various techniques to regularize the delta function [24,25]. In this tutorial we choose:
where is the regularization parameter that controls the support of the regularized delta functions . Typically should have a size comparable to the mesh spacing .
Set up the regularized delta functions centered at .
In[55]:=
Click for copyable input
Determine the appropriate mesh spacing using the function introduced in Model Parameter Setup.
In[56]:=
Click for copyable input
Set the monopole source strength and the monopole source . Here we define the regularization parameter to be half of the mesh spacing: .
In[57]:=
Click for copyable input
Set up a 1D acoustic model with the monopole source term centered at .
In[59]:=
Click for copyable input
Apply absorbing boundary conditions on both sides of the domain to avoid reflection. Then the PDE for the sound pressure field is given by:
In[61]:=
Click for copyable input
Solve the PDE with ParametricNDSolveValue.
In[62]:=
Click for copyable input
As shown in the section entitled Time Harmonic Analysis, the solution from the Helmholtz PDE is the complex valued amplitude function , which contains the information of both phase and amplitude. The amplitude of the sound pressure corresponds to the absolute value or complex modulus .
Inspect the amplitude of the sound pressure field at frequencies
In[63]:=
Click for copyable input
Out[65]=
In a 1D domain the sound pressure amplitude is constant for the monopole source at each particular frequency. Note that for a given monopole source strength , the pressure amplitude reduces with higher frequencies. The analytical solution [26] for the amplitude in 1-dimension is given by
It is more intuitive to consider the sound pressure as a function of time. To do so the amplitude function should be transformed via the harmonic wave relation (27): .
Inspect the monopole source in the time domain and compare with the analytical amplitudes.
In[66]:=
Click for copyable input
Out[67]=
The monopole source is a point source that radiates sound isotropically at . The blue line here is the transient sound pressure , and the gray line is the analytical pressure amplitude inserted for a visual verification.
As a comparison, a 2D monopole source is constructed in a similar manner.
Define a 2D simulation domain.
In[68]:=
Click for copyable input
Set up a 2D acoustic model with the monopole source term centered at .
In[70]:=
Click for copyable input
Apply an absorbing boundary condition on the outer boundary to avoid reflection, then the PDE for the sound pressure field is given by:
In[71]:=
Click for copyable input
Solve the PDE with ParametricNDSolveValue.
In[73]:=
Click for copyable input
Set up a legend bar and ContourPlot options for scaling from to .
In[74]:=
Click for copyable input
Inspect the amplitude of the sound pressure field at frequency .
In[76]:=
Click for copyable input
Out[77]=
For 2D and 3D monopoles, as the radiated sound wave spreads out from the source a wider wavefront will be formed. Therefore, with a given monopole strength the pressure amplitude will decrease with the distance to the source location

Dipole Sources

A dipole source consists of two monopole sources of equal strength but opposite phase and separated by a distance . The separation distance is small compared to the wavelength of sound . An example of an acoustic dipole would be a small rigid sphere that oscillates sinusoidally [28]. Unlike a single monopole source , a dipole source does not radiate sound isotropically.
In the acoustics time domain tutorial it has been shown that a dipole source could be modeled by two identical monopole sources . Here we directly jump to model a dipole source with the source term .
To make use of a dipole source we have to specify the dipole source strength and the source location . The dipole source term may be written as
where is the regularize delta function at the source location.
Set the dipole source strength and the 1D dipole source . Here we define the regularization parameter to be half of the mesh spacing: .
In[78]:=
Click for copyable input
Set up a 1D acoustic model with a dipole source term centered at .
In[80]:=
Click for copyable input
With absorbing boundary conditions on both sides the PDE for the sound pressure field is given by:
In[82]:=
Click for copyable input
Solve the PDE with ParametricNDSolveValue.
In[83]:=
Click for copyable input
As shown in the section entitled Time Harmonic Analysis, the solution from the Helmholtz PDE is the complex valued amplitude function , which contains the information of both phase and amplitude. The amplitude of the sound pressure corresponds to the absolute value or complex modulus .
Inspect the amplitude of the sound pressure field at frequencies
In[84]:=
Click for copyable input
Out[86]=
In 1D domain the sound pressure amplitude is constant outside each dipole source region. Since the embedded monopole sources are in opposite phase, the amplitude at the source location sums up to zero. The analytical solution [29] for the pressure amplitude in 1D is given by
The numerical solution plot above, however, is slightly different from the analytical solution due to the discrete nature of the regularized delta function used. A finer mesh can be used to reduce the numerical error.
It is more intuitive to consider the sound pressure as a function of time. To do so the amplitude function should be transformed via the harmonic wave relation (30): .
Inspect the dipole source in the time domain and compare with the analytical amplitudes.
In[87]:=
Click for copyable input
Out[88]=
The dipole source does not radiate sound isotropically. The resultant sound waves are sinusoidal but in opposite phase.
As a comparison, a 2D dipole source is constructed in a similar manner.
Set the 2D dipole source with a directivity angle .
In[89]:=
Click for copyable input
Set up a 2D acoustic model with the dipole source term centered at .
In[91]:=
Click for copyable input
Apply an absorbing boundary condition on the outer boundary to avoid reflection, then the PDE for the sound pressure field is given by:
In[94]:=
Click for copyable input
Solve the PDE with ParametricNDSolveValue.
In[96]:=
Click for copyable input
Set up a legend bar and ContourPlot options for scaling from to .
In[97]:=
Click for copyable input
Inspect the amplitude of the sound pressure field at frequency .
In[99]:=
Click for copyable input
Out[100]=
For 2D and 3D dipoles the radiated sound wave does not spread out equally in all directions [31]. Therefore, with a given dipole strength the pressure amplitude will depend on both the spatial direction and the distance to the source location.
Sound Propagation in Lossy Media
As shown in the acoustics time domain tutorial, it is possible to model sound propagation in lossy media with a given porosity , flow resistivity and an effective bulk modulus . The modified wave equation is
Inserting the harmonic wave relation and factoring out the common term , the equation (32) becomes
Equation (33) is the modified Helmholtz equation that is used to model sound attenuation in the frequency domain.
The effective bulk modulus is frequency-dependent [34], and can be approximated by an empirical formula:
To illustrate the lossy media model we consider a sound wave propagating through a tube filled with a porous absorber. The porosity and flow resistivity of this material are given by and , respectively.
Set up a 1D source free acoustic model with an attenuation term.
In[101]:=
Click for copyable input
Calculate the effective bulk modulus .
In[102]:=
Click for copyable input
An absorbing boundary condition is added on the right to avoid a reflection of the wave. Note that the NeumannValue is set to [35] to accommodate for the modified Helmholtz equation (36).
Set up an absorbing boundary condition for the lossy media model.
In[103]:=
Click for copyable input
A harmonic sound pressure source with amplitude is added on the left end using DirichletCondition.
In[104]:=
Click for copyable input
The PDE for the sound pressure field is given by:
In[106]:=
Click for copyable input
Solve the PDE with ParametricNDSolveValue.
In[107]:=
Click for copyable input
Inspect the exponential decay of the pressure amplitude at frequencies with a logarithmically scaled axis.
In[108]:=
Click for copyable input
Out[110]=
Since the effective bulk modulus in (37) is a frequency-dependent variable, the attenuation rate of the amplitude also varies with different frequencies. The small wiggles near the right end result from the numerical reflection on the absorbing boundary, and can be reduced by using a finer mesh.
Inspect the sound wave in the time domain.
In[111]:=
Click for copyable input
Out[111]=PlayAnimation
The amplitude of the sound pressure wave decays as it propagates to the right. Note that the attenuation rate slightly varies with different frequencies.
A Comparison of Time Domain and Frequency Domain Modeling
An acoustic system can be modeled both in the time domain and the frequency domain. As shown in the acoustics time domain tutorial, the wave equation is used to find a transient solution of a sound wave in the time domain. Whenever the excitation of a sound pressure field is time-harmonic, the Helmholtz equation can be used to directly solve for a steady-state solution of a sound pressure field in the frequency domain. Each approach has its own strength and constraints, and is summarized in the following table.
    

Time Domain Modeling Frequency Domain Modeling
Governing PDEWave equationHelmholtz equation
Dependent variabletransient p(t,X)stationary p(X)
Harmonic excitationyesyes
Inharmonic excitationyesno
Computational costhighlow
Accuracylowhigh

In other words, if an acoustic system has a time inharmonic dependency then it needs to be modeled in the time domain, and the benefits of using the Helmholtz equation can not be exploited. But in all other cases making use of the model in the frequency domain is beneficial.
To illustrate this behaviour an acoustic system of a right-traveling wave is considered in the next example. The simulation is done once in the time domain, and once in the frequency domain. In both models the simulation domain is set as four times of the wavelength , and the frequency is arbitrarily chosen at .
Set the sampled frequency , angular frequency and the wavelength .
In[112]:=
Click for copyable input
A radiation boundary condition is added on the left to produce a harmonic sound wave, and an absorbing boundary condition is placed on the right to avoid wave reflection.
Set the boundary conditions for time-domain and frequency-domain models.
In[115]:=
Click for copyable input
In[118]:=
Click for copyable input
Next we will analyse the model in the time domain, and following that a frequency domain analysis will be performed. The results will be compared to each other.

Wave Equation: Time Domain Modeling

For the wave equation model, the simulation end time is defined as the time required for the sound pressure field to reach its dynamic steady state.
Set the period and the simulation end time .
In[121]:=
Click for copyable input
Set up a 1D wave equation for the time domain modeling.
In[123]:=
Click for copyable input
Out[123]=
The initial conditions are set as an undisturbed domain.
In[124]:=
Click for copyable input
The PDE for the sound pressure field is given by:
In[125]:=
Click for copyable input
Solve the wave equation and monitor time/memory usage.
In[126]:=
Click for copyable input
Inspect the sound wave in the time domain
In[128]:=
Click for copyable input
Out[130]=
The result shows a transient sound pressure wave traveling to the right.
Next, we build the same model with the Helmholtz equation.

Helmholtz Equation: Frequency Domain Modeling

Inspect the setting of a 1D Helmholtz equation for the frequency domain modeling.
In[131]:=
Click for copyable input
Out[131]=
The PDE for the sound pressure field is given by:
In[132]:=
Click for copyable input
Solve the Helmholtz PDE and monitor time/memory usage.
In[133]:=
Click for copyable input
Since the Helmholtz equation directly computes a stationary solution of a sound pressure field, there is no need to do the time integration in the solving process. The computational cost is thus reduced significantly.
Inspect the steady-state solution of a sound wave.
In[135]:=
Click for copyable input
Out[135]=
The result of the Helmholtz equation is a steady-state sound pressure field . To transform the solution into the time domain the harmonic wave relation can be used:
To transform a transient solution into the frequency domain, however, requires a Fourier Analysis [38]. The process is to decompose the transient sound signal into a sum of harmonic signals. Each harmonic signal has a specific frequency and a relative magnitude, which allows us to map the transient signal into the frequency domain.
The computational cost for solving the Helmholtz equation is so much lower that it is possible to solve the Helmholtz PDE repetitively with different frequencies, which makes it suitable for frequency domain modeling.
Next, we compare the accuracy of the wave equation model and the Helmholtz equation model.

Accuracy Comparison

Set up the analytical solution.
In[136]:=
Click for copyable input
Compare the error of the wave equation model and the Helmholtz equation model.
In[137]:=
Click for copyable input
Out[137]=
Since the wave equation model is a time-dependent PDE while the Helmholtz equation model is not, the former is subject to an extra error from the numerical time integration. The Helmholtz model is therefore more accurate than the wave equation model.
Acoustic Boundary Conditions
In the acoustics time domain tutorial, we have shown the details of common acoustic boundary conditions and how they can be modeled in the time domain. To avoid the repetition the following section only describes how to build these boundary conditions in the frequency domain. For readers who are interested in the derivation and the physical explanation, please refer to the acoustics time domain tutorial.
Most common boundary conditions in acoustics can be modeled with DirichletCondition and NeumannValue, and can be categorized in the following three types:
For each boundary condition we will state weather the particular boundary condition is applicable for Time Harmonic Analysis or Eigenfrequency Analysis or both in the following manner:

Analysis Type

Applicable

Yes/No

Yes/No

Impedance Boundary Conditions
When a sound wave transits to another medium or encounters a partially reflective boundary, for example a porous surface, part of the sound wave will be reflected at the interface and part will be transmitted across the interface.
One property that is used to formulate the relation between the incident, reflected and transmitted wave is called the specific impedance. The specific impedance is the ratio of the sound pressure to the sound particle velocity as the wave moves through the medium, which is defined as:
From the acoustics time domain tutorial we know that the impedance boundary condition can be formulated with a given boundary impedance as
Inserting the harmonic wave relation and factoring out the common term , then the impedance boundary condition can be applied in the frequency domain as
Here the wave number denotes the ratio of angular frequency to the speed of sound.
An impedance boundary condition can be used with:

Analysis Type

Applicable

Yes

No

Impedance Boundary Conditions in Time Harmonic Analysis

To illustrate the impedance boundary condition a tube with a porous surface at the right end is considered in the next example. The porous surface is treated as an impedance boundary since it is a partially reflective boundary.
The boundary impedance of the porous surface is assumed to be for demonstration purposes.
In[138]:=
Click for copyable input
Set up a impedance boundary condition on the right end with NeumannValue.
In[139]:=
Click for copyable input
With a radiation boundary condition on the left end the PDE for the sound pressure field is given by:
In[140]:=
Click for copyable input
Solve the PDE with ParametricNDSolveValue.
In[141]:=
Click for copyable input
Inspect the amplitude of the sound pressure field at frequencies
In[142]:=
Click for copyable input
Out[144]=
It is more intuitive to inspect the sound wave in the time domain. Recall that the solution can be transformed back into time domain with the harmonic wave relation (39): .
Inspect the sound wave in time domain
In[145]:=
Click for copyable input
Out[145]=PlayAnimation
Since the impedance boundary models a partially reflective boundary, there is more energy moving to the right than there is being reflected, which makes the resulting wave appear to travel to the right. Note that the maximum and the minimum values of the pressure amplitude are fixed in space (dashed lines). This type of wave is called a partial standing wave.
The ratio between the maximum and the minimum amplitudes is known as the standing wave ratio (SWR).
Calculate the standing wave ratio (SWR) at frequencies
In[146]:=
Click for copyable input
Out[148]=
The standing wave ratio (SWR) is shown to be independent of the frequency. For an unknown boundary the SWR can be used to measure the impedance and the reflection coefficient given the specific impedance, , of the domain.
Here and are the amplitude of the reflected and the incident wave, respectively.
Absorbing Boundary Conditions
Typically a simulation domain that extends to infinity is not a feasible option in a simulation. Absorbing boundary conditions are a methodology used to model infinite domains. An absorbing boundary condition (ABC) works by absorbing an incoming wave and thus makes the model behave as if it had infinite extent. ABC are not the only way to model simulation domains with infinite extent. Perfectly matched layers (PML) may be used as an alternative approach to an ABC.
From the time domain tutorial we know that the absorbing boundary condition is given by
and inserting the harmonic wave relation gives the absorbing boundary condition in the frequency domain.
An absorbing boundary condition can be used with:

Analysis Type

Applicable

Yes

No

Absorbing Boundary Conditions in Time Harmonic Analysis

As an example we look at an infinitely long tube with the computational domain set from to . To model the continuation of the domain we add a plane wave absorbing boundary condition at the right end.
Inspect the setting of an absorbing boundary condition on the right end.
In[149]:=
Click for copyable input
Out[149]=
With a radiation boundary condition on the left end the PDE for the sound pressure field is given by:
In[150]:=
Click for copyable input
Solve the PDE with ParametricNDSolveValue.
In[151]:=
Click for copyable input
Inspect the amplitude of the sound pressure field at frequencies
In[152]:=
Click for copyable input
Out[154]=
Since the resulting sound pressure field is simply a right-traveling harmonic wave, the pressure amplitude is fixed at throughout the domain for all frequencies.
Inspect the sound wave in the time domain.
In[155]:=
Click for copyable input
Out[155]=PlayAnimation
The incoming wave is absorbed at the right hand boundary as if the simulation domain had infinite extent.
Sound Hard Boundary Conditions - Walls
On a sound hard boundary, the normal component of the sound particle velocity is zero since no forward motion is possible.
Substituting (40) into the momentum conservation equation and applying the harmonic wave relation (41) , then the sound hard boundary condition can be formulated in the frequency domain as
Set up a sound hard boundary condition on the right end with NeumannValue.
In[156]:=
Click for copyable input
If no boundary condition is specified on any part of the boundary then by default a Neumann zero boundary condition is implicitly used. This implies that a sound hard boundary is the default boundary condition used if no boundary condition is specified at a given boundary.
In[157]:=
Click for copyable input
A sound hard boundary condition can be used with:

Analysis Type

Applicable

Yes

Yes

Sound Hard Boundary Conditions in Time Harmonic Analysis

As an example for a time harmonic analysis, we look at a tube with one end closed.
With a radiation boundary condition on the left end the PDE for the sound pressure field is given by:
In[158]:=
Click for copyable input
Solve the PDE with ParametricNDSolveValue.
In[159]:=
Click for copyable input
Inspect the amplitude of the sound pressure field at frequencies
In[160]:=
Click for copyable input
Out[162]=
The shape of the amplitude field shows several minima and maxima points, which are known as "nodes" and "antinodes", respectively. In the time domain nodes are positions where the standing wave has no displacement and antinodes are the positions with maximal displacement. In the frequency domain the displacement manifest themselves as minima and maxima. Since no forward motion is possible on the sound hard boundary, the sound pressure is fixed at its maximum at the right end, which means it is one of the antinodes. The maximum value corresponds to the double of the amplitude set by the radiation boundary.
It is more intuitive to inspect the result in the time domain.
Inspect the sound wave in the time domain.
In[163]:=
Click for copyable input
Out[163]=PlayAnimation
Note that the sound wave neither moves right nor left but simply oscillates in time. This type of wave is known as a standing wave, which is formed by the superposition of two waves traveling in opposite directions.
In this case the right-traveling wave is produced by the radiation boundary and the left-traveling wave is the reflected wave from the sound hard boundary. For readers who are interested in the way that travelling waves superimpose to give a standing wave, please refer to the acoustics time domain tutorial.

Sound Hard Boundary Conditions in Eigenfrequency Analysis

Unlike a time harmonic analysis, an eigenfrequency analysis aims to find the eigenfrequencies and the corresponding eigenmodes (eigenfunctions) of an acoustic system. In the next example we consider a tube with both end closed.
Set up sound hard boundary conditions on both ends.
In[164]:=
Click for copyable input
Solve for the five smallest eigenvalues and eigenmodes with NDEigensystem.
In[166]:=
Click for copyable input
Calculate the corresponding eigenfrequencies with the relation .
In[167]:=
Click for copyable input
Out[167]=
Note that the first eigenfrequency is zero, corresponding to the solution without any sound. The first pair of eigenfrequency/eigenmode is therefore a trivial solution, and is denoted as eigenmode 0.
Inspect the eigenmode within the domain.
In[168]:=
Click for copyable input
Out[168]=PlayAnimation
As explained previously, the sound pressure is fixed at its maximum on the sound hard boundaries since no forward motion is possible.
Normal Velocity Boundary Conditions
When a nonzero, time harmonic sound particle velocity is specified at a boundary, then this type of boundary is called a normal velocity boundary.
Substituting (42) into the momentum conservation equation and applying the harmonic wave relation (43) , then the normal velocity boundary condition can be formulated in the frequency domain as
A normal velocity boundary condition can be used with:

Analysis Type

Applicable

Yes

No

Normal Velocity Boundary Conditions in Time Harmonic Analysis

In the following example a harmonic vibration is introduced at the right hand boundary and vibrates with a known velocity amplitude . The sound field can be calculated by using NeumannValue as shown below.
Set up a normal velocity boundary condition on the right end with NeumannValue.
In[169]:=
Click for copyable input
With an absorbing boundary condition on the left end the PDE for the sound pressure field is given by:
In[171]:=
Click for copyable input
Solve the PDE with ParametricNDSolveValue.
In[172]:=
Click for copyable input
Inspect the amplitude of the sound pressure field at frequencies
In[173]:=
Click for copyable input
Out[175]=
Since the resulting wave is simply a left-traveling harmonic wave, the amplitude distribution is fixed throughout the domain for all frequencies.
Inspect the sound wave in the time domain.
In[176]:=
Click for copyable input
Out[176]=PlayAnimation
The normal velocity boundary generates a harmonic wave on the right end that propagates to the left.
As shown in the acoustics time domain tutorial, a normal velocity boundary could be replaced by a pressure source boundary when the specific impedance is known.
Sound Soft Boundary Conditions
On a sound soft boundary the medium pressure is set equal to an ambient reference pressure, which means that the sound pressure at the boundary is fixed at zero: .
A sound soft boundary condition can be used with:

Analysis Type

Applicable

Yes

Yes

Sound Soft Boundary Conditions in Time Harmonic Analysis

As an example for the time harmonic analysis, we look at a tube with one end open. The open-ended side is treated as a sound soft boundary since there is no constraint to limit the sound wave movement.
Set up a sound soft boundary condition on the right end with DirichletCondition.
In[177]:=
Click for copyable input
With a radiation boundary condition on the left end the PDE for the sound pressure field is given by:
In[178]:=
Click for copyable input
Solve the PDE with ParametricNDSolveValue.
In[179]:=
Click for copyable input
Inspect the amplitude of the sound pressure field at frequencies
In[180]:=
Click for copyable input
Out[182]=
With the sound pressure is a minimum at , where a sound soft boundary is positioned, and is thus called a node. The maximum value of the sound pressure field corresponds to the double of the amplitude set by the radiation boundary.
Inspect the sound wave in the time domain.
In[183]:=
Click for copyable input
Out[183]=PlayAnimation
Similar as the sound hard boundary condition, the resulting wave is a standing wave that is formed by the superposition of two waves traveling in opposite directions.
In this case the right-traveling wave is produced by the radiation boundary and the left-traveling wave is the reflected wave from the sound soft boundary. For readers who are interested in the way that travelling waves superimpose to give a standing wave, please refer to the acoustics time domain tutorial.

Sound Soft Boundary Conditions in Eigenfrequency Analysis

To illustrate the behavior of sound soft boundaries in an eigenfrequency analysis, we consider a tube with both ends open.
Set up sound soft boundary conditions on both ends.
In[184]:=
Click for copyable input
Solve the five smallest eigenvalues and eigenmodes with NDEigensystem.
In[186]:=
Click for copyable input
Calculate the corresponding eigenfrequencies with the relation .
In[187]:=
Click for copyable input
Out[187]=
Inspect the eigenmodes within the domain.
In[188]:=
Click for copyable input
Out[188]=PlayAnimation
As expected, the sound pressure on the sound soft boundaries is fixed at zero.
Pressure Source Boundary Conditions
We speak of a pressure source boundary condition when a nonzero, time harmonic sound pressure is specified at a boundary. Both DirichletCondition and NeumannValue can be used to specify a pressure source boundary condition.
A pressure source boundary condition can be used with:

Analysis Type

Applicable

Yes

No

Pressure Source Boundary Conditions in Time Harmonic Analysis

Dirichlet Model
To formulate the Dirichlet condition of a pressure source boundary, we simply rewrite equation (44) with the harmonic wave relation (45) :
Here denotes the amplitude of the pressure source.
Set up a pressure source boundary condition on the right end with DirichletCondition.
In[189]:=
Click for copyable input
With an absorbing boundary condition on the left end the PDE for the sound pressure field is given by:
In[191]:=
Click for copyable input
Solve the PDE with ParametricNDSolveValue.
In[192]:=
Click for copyable input
Inspect the amplitude of the sound pressure field at frequencies
In[193]:=
Click for copyable input
Out[195]=
Since the resulting wave is simply a left-traveling harmonic wave, the amplitude distribution is fixed at throughout the domain.
Inspect the sound wave in the time domain.
In[196]:=
Click for copyable input
Out[196]=PlayAnimation
Similar to the normal velocity boundary, the pressure source generates a harmonic wave on the right end that propagates to the left.
Neumann Model
A pressure source can also be modeled with a NeumannValue. As shown in the acoustics time domain tutorial the NeumannValue setting for a pressure source is given by
Insert the harmonic wave relation (46) , then the pressure source boundary condition can be formulated in the frequency domain as:
Set up a pressure source boundary condition on the right end with NeumannValue.
In[197]:=
Click for copyable input
With an absorbing boundary condition on the left end the PDE for the sound pressure field is given by:
In[199]:=
Click for copyable input
Solve the PDE with ParametricNDSolveValue.
In[200]:=
Click for copyable input
Inspect the amplitude of the sound pressure field at frequencies
In[201]:=
Click for copyable input
Out[203]=
Inspect the sound wave in the time domain.
In[204]:=
Click for copyable input
Out[204]=PlayAnimation
The animation shows the same effect as the Dirichlet model. For readers who are interested in the trade-off between the Neumann model and the Dirichlet model of a pressure source boundary, please refer to the corresponding section in the acoustics time domain tutorial.
Radiation Boundary Conditions
A radiation boundary is a hybrid boundary condition that combines the properties of a pressure source boundary and an absorbing boundary. A time harmonic sound pressure is specified at a boundary to produce an incoming wave, however, unlike a simple pressure source the radiation boundary allows an outgoing wave to leave the computational domain with little reflection.
As shown in the acoustics time domain tutorial the radiation boundary condition is given by
Insert the harmonic wave relation (47) , then the radiation boundary condition can be formulated in the frequency domain as:
A radiation boundary condition can be used with:

Analysis Type

Applicable

Yes

No

Radiation Boundary Conditions in Time Harmonic Analysis

As an example a semi-infinite tube with one end closed is considered in the next example. The right end is treated as a radiation boundary where the harmonic sound wave enters the computational domain, and the closed left end is set as a sound hard boundary.
Inspect the setting of a radiation boundary condition on the right end. By default a sound hard boundary condition is implicitly used at the left end.
In[205]:=
Click for copyable input
Out[205]=
The PDE for the sound pressure field is given by:
In[206]:=
Click for copyable input
Solve the PDE with ParametricNDSolveValue.
In[207]:=
Click for copyable input
Inspect the amplitude of the sound pressure field at frequencies
In[208]:=
Click for copyable input
Out[210]=
Inspect the sound wave in the time domain.
In[211]:=
Click for copyable input
Out[211]=PlayAnimation
The resulting sound pressure field shows a standing wave, which is formed by the superposition of both incoming and outgoing waves. In this case the radiation boundary at generates an incoming left-traveling wave, and the implicit sound hard boundary at reflects the wave to the right as an outgoing wave, which leaves the domain from the radiation boundary without constrained. For readers who are interested in the way that travelling waves superimpose to give a standing wave, please refer to the acoustics time domain tutorial.
Perfectly Matched Layer (PML)
A perfectly matched layer (PML) is a method to model simulation domains with infinite extent. As such a PML is an alternative method to absorbing boundary conditions . The following section shows how to implement a PML for a Helmholtz PDE used in frequency domain modeling. The full derivation and the explanation can be found in the appendix section of the acoustics time domain tutorial.
To implement a PML two things need to happen. First the simulation domain needs to be enlarged. This extension is the region in which the PML is active. Second, a coordinate transformation of the PDE is done.
The 3 dimensional Helmholtz equation after the PML coordinate transformation is given in by [48]
Here , and are the absorbing coefficients of the PML. Three auxiliary variables , and are introduced to control the PML attenuation in each dimension.
In the 1D case where , the equation (49) simplifies to
The absorption coefficient is a tuning parameter to be chosen, and is set to increase linearly within the PML layer from to . A thinner PML requires a greater value of but large tends to increase numerical reflections.
The perfectly matched layer can be applied in Time Harmonic Analysis only.
Perfectly Matched Layer in Time Harmonic Analysis
As an example we consider a simulation of a 1D domain from to using a computational domain that ranges from to with a PML at the right from to .
Define the 1D domain and the PML width.
In[212]:=
Click for copyable input
Define the PML model with an auxiliary variables .
In[214]:=
Click for copyable input
Out[214]=
Specify and set the absorption coefficient to increase linearly within the PML layer.
In[215]:=
Click for copyable input
Inspect the constructed absorption coefficient .
In[217]:=
Click for copyable input
Out[217]=
As seen in the original non-PML region.
With a radiation boundary condition on the left end and a specified auxiliary variables , then the Helmholtz PDE is given by:
In[218]:=
Click for copyable input
Solve the PDE with ParametricNDSolveValue.
In[219]:=
Click for copyable input
Visualize the amplitude of the sound pressure field including the PML region. The sampled frequency is taken as .
In[220]:=
Click for copyable input
Out[222]=
The sound pressure in the non-PML region is fixed at the incident amplitude , and then decays to zero within the PML region. To minimize the numerical reflection at the PML interface, the absorption coefficient and the PML width should be chosen carefully. The tradeoff between these two parameters is discussed in the acoustics time domain tutorial.
Inspect the wave propagation in the time domain.
In[223]:=
Click for copyable input
Out[223]=PlayAnimation
The wave is attenuated within the PML but remains unchanged in the original domain, and the attenuation rate is independent of frequencies.
Nomenclature
    

SymbolDescriptionUnit
ρdensity of a medium[kg/m3]
cspeed of sound in a medium[m/s]
psound pressure[Pa]
plocal sound amplitude[Pa]
specified boundary pressure[Pa]
conjugate of sound pressure[Pa]
ttime[s]
tendsimulation end time[s]
Xposition vector[m]
sdirection switchN/A
Foptional dipole source[N/m3]
dipole source strength[N/m3]
θdipole directivity angle[rad]
Qoptional monopole source[1/s2]
monopole source strength[1/s2]
Nullseperation distance of dipole source[m]
λwavelength of sound[m]
Ω simulation domain[m]
Nullwave number[rad/m]
Nullsound wave frequency[Hz]
ωsound wave angular frequency[rad/s]
δDirac delta functionN/A
regularized delta functionN/A
Nullregularization parameter[m]
Nullmesh spacing[m]
Nullsound source location[m]
Nulleffective bulk modulus[Pa]
αattenuation factor[m2/(s·N)]
ϕporosityN/A
Vvviod volume[m3]
VTtotal volume[m3]
Rfflow resistivity[kg/(m3·s)]
βstandard deviation of a Gussian pulse[m]
ζsound particle displacement[m]
Nullsound particle velocity[m/s]
Nullspecified boundary velocity[m/s]
Tsound wave period[t]
Zcharacteristic impedence[Null]
Zbboundary impedance[Null]
Aramplitude of reflected wave[Pa]
Aiamplitude of incident wave[Pa]
σabsorbtion coefficient of PML[rad/(s·m)]
σmaxmaximum value of absorbtion coefficeint[rad/(s·m)]

References

1.  Ihlenburg, Frank. The Medium-Frequency Range in Computational Acoustics: Practical and Numerical Aspects. Journal of Computational Acoustics, Vol.11, No. 2 175-193, 2003.

2.  Heutschi, Kurt. Lecture Notes on Acoustics I. Swiss Federal Institute of Technology Zurich, 2016.

3.  Johnson, Steven. Notes on Perfectly Matched Layers (PMLs). MIT, 2010.

4.  Bilbao, Stefan and Hamilton, Brian. Directional Source Modeling In Wave-Based Room Acoustics Simulation. IEEE, 2017.

5.  Peskin, Charles. The Immersed Boundary Method. Cambridge University, 2002.

6.  Russell, Daniel, Titlow, Joseph and Bemmen, Ya-Juan. Acoustic monopoles, dipoles and quadropoles: An experiment revisited. American Journal of Physics 67, 660, 1999.

7.  Vita, Micro. The Wave Equation with a Source. Oklahoma State University.

8.  J. De Moerloose and M. A. Stuchly, Behavior of Berenger's ABC for evanescent waves, IEEE Microwave and Guided Wave Letters, vol. 5, no. 10, pp. 344-346, Oct. 1995.

9.  J. Berenger, Evanescent waves in PML's: origin of the numerical reflection in wave-structure interaction problems, IEEE Transactions on Antennas and Propagation, vol. 47, no. 10, pp. 1497-1503, Oct. 1999.

10.  J. De Moerloose, Jan & Stuchly, Maria. Reflection analysis of PML ABCs for low-frequency applications, IEEE Microwave and Guided Wave Letters, vol. 6., no. 4, pp. 177-179, Apr. 1996.

11.  E. Turkel and A. Yefet. Absorbing PML boundary layers for wave-like equations, Applied Numerical Mathematics, vol. 27, pp. 533-557, 1998.

12.  G. Pan, A. Abubakar and T. Habashy. An effective perfectly matched layer design for acoustic fourth-order frequency-domain finite-difference scheme, Geophysical Journal International, vol. 188, pp. 211-222, 2012.

13.  T. J. Cox, P. D'Antonio. Acoustic absorbers and diffusers: Theory, design, and application, London: Spon Press, 2004.

14.  E. W. Weisstein. Fast Fourier Transform, MathWorld--A Wolfram Web Resource. Retrieve from: http://mathworld.wolfram.com/FastFourierTransform.html.