Acoustics in the Time Domain

Introduction
Acoustics is the field of physics that models sound by changes in pressure. The changes in pressure are described by a wave equation. This tutorial gives an introduction to modeling sound with the wave equation in the time domain and presents various aspects of the modeling process. The modeling process results in partial differential equation (PDE) models that are solved with NDSolve. Furthermore, in this tutorial different types of sound sources are introduced as well as an overview of how various real world sound barriers can be modeled with the available PDE boundary conditions.
In this tutorial time domain modeling is used to investigate how sound waves behave in time. It should be noted that sound can also be modeled in the frequency domain. Employing the frequency domain for modeling illustrates how sound waves are distributed over a range of frequencies. Frequency domain modeling is widely used when analyzing acoustic systems with a frequency dependency, such as noise filters.
Details about the modeling of sound in the frequency domain can be found in a separate tutorial entitled Acoustics in the Frequency Domain.
The symbols and corresponding units used throughout this tutorial are summarized in the Nomenclature section.
Load the finite element package.
In[23]:=
Click for copyable input
Wave Equation
Introduction to the Wave Equation
The wave equation (1) is used for modeling the time dependent (transient) propagation of sound in a lossless medium.
The dependent variable in the wave equation is the sound pressure . The sound pressure is propagating in a medium with density at the speed of sound . Sound pressure can be understood as the local pressure deviation from an ambient reference pressure: . The pressure varies with time and position . Terms and are monopole and dipole sources, respectively. The Source Types section describes these sound sources.
Sound can propagate through all forms of matter without a net transport of mass and is typically transmitted as longitudinal waves by alternating the sound pressure of the medium along the direction of wave propagation. For solid media sound will be transmitted longitudinally and additionally transversely such that waves also oscillate perpendicular to the direction of propagation [2]. The following animation illustrates particle motions of longitudinal waves and transverse waves [3].

The animation illustrates how the highlighted particle oscillates horizontally on longitudinal waves, but fluctuates vertically on transverse waves.
Defining a function that represents the spatial terms of an acoustic model will make the setup for the wave equation more convenient.
Define the acoustic model representation function.
In[25]:=
Click for copyable input
Set up a 1D time dependent source free acoustic model.
In[28]:=
Click for copyable input
Out[28]=
Various sections in the documentation explain the use of inactive PDE operators. Please refer to "Numerical Solution of Partial Differential Equations".
To make use of specific material parameters in the equation we extract relevant data from the ThermodynamicData.
In[29]:=
Click for copyable input
The actual numerical material parameters can be put into the model using replacement:
In[31]:=
Click for copyable input
Out[31]=
Model Parameter Setup
Sound waves can be classified into two categories based on their form of time dependence. Time harmonic waves have sinusoidal time variation with a single frequency. Time inharmonic waves have arbitrary time variation with multiple frequencies.
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 max cell measure for both time harmonic and time inharmonic waves.
Throughout this tutorial we will use a particular example frequency for time harmonic waves and a Gaussian pulse as an example of time inharmonic waves.
The default resolution of the mesh is set to 20 elements per wavelength , which means the max cell measure would be twenty times smaller than the sound wavelength . Typically this is sufficient to resolve the waves accurately [4]. 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[32]:=
Click for copyable input
The wavelength is computed differently for time harmonic and time inharmonic waves. For time harmonic waves or , where the frequency is related to the angular frequency by . And for time inharmonic waves . Here, is the standard deviation of the Gaussian pulse [5].
Define a function named AcousticOptions that gives NDSolve options appropriate for solving acoustic problems.
In[33]:=
Click for copyable input
The following model parameters are used for the examples in this tutorial. These parameters define the simulation domain and the simulation end time .
Set up model parameters, including the angular frequency for time harmonic examples.
In[35]:=
Click for copyable input
If an acoustic wave begins with an undisturbed domain, the initial sound pressure and its derivative are set to zero.
Define
Click for copyable input
to represent an undisturbed initial condition.
In[40]:=
Click for copyable input
In the examples of time inharmonic waves we will typically be using Gaussian pulses with standard deviation .
Define
Click for copyable input
to represent a Gaussian pulse:
In[41]:=
Click for copyable input
In many examples we will be using an absorbing boundary condition (ABC) to avoid wave reflection.
Define
Click for copyable input
and
Click for copyable input
to be absorbing boundary conditions on the left end and
right end, respectively.
In[44]:=
Click for copyable input
Wave Types
Sound waves can be classified into two categories based on their form of time dependence. Time harmonic waves have sinusoidal time variation with a single frequency. Time inharmonic waves have arbitrary time variation with multiple frequencies.

Time Harmonic Waves

A wave depending on a single frequency and a sinusoidal variation in time is denoted as a time harmonic wave and is given by (6).
Here represents a spatial distribution of the sound pressure amplitude.
An arbitrary wave can be decomposed into an infinite sum of harmonic waves with discrete frequencies [7]. This process is called frequency analysis. Most acoustical analysis involves frequency analysis since many acoustic systems, such as human ears or noise controllers, are frequency dependent. Harmonic wave behaviour, the basis of frequency analysis, is a fundamental topic in acoustical engineering.
Time harmonic waves are often expressed in a complex exponential representation:
For a wave with harmonic space dependency the amplitude function is:
where is the amplitude of the harmonically oscillating quantity and is the wave number and is interpreted as phase change per unit distance.
The variable in the exponent of defines the direction of travel of the wave. Setting makes the wave propagate in the negative direction and in the positive x direction.
Combining (8) and (9) a traveling time harmonic wave is expressed as:
Use Manipulate to investigate a wave with adjustable frequency and wave number . Assume that the amplitude .
In[46]:=
Click for copyable input
Out[46]=PlayAnimation
There are four model parameters in the animation; a time slider, a direction switch , a frequency and wave number selector. The number of phase changes per unit time increase when selecting a higher frequency . The number of phase changes per unit distance increase as the wave number increase. The wave propagation can be observed by playing the time slider.
When is substituted in the source free wave equation (10), the equation simplifies to an inhomogeneous Helmholtz equation that is a second order time independent differential equation.
In the case where there are no sources this PDE system becomes an eigenvalue problem that can be solved with NDEigensystem. The set of eigenvalues then determines the natural frequencies of the corresponding acoustic system.
The details of the sound modeling with the Helmholtz equation can be found in a separate tutorial entitled Acoustics in the Frequency Domain.
Next, we build a 1D time harmonic wave model using the wave equation (11) and solve it numerically using NDSolveValue.
A sinusoidal source on the left end of will produce a time harmonic wave. The source can be added using DirichletCondition.
In[47]:=
Click for copyable input
The PDE for the sound pressure field is given by:
In[48]:=
Click for copyable input
Substitute the actual material parameters in the PDE.
In[49]:=
Click for copyable input
Solve the PDE with NDSolveValue.
In[50]:=
Click for copyable input
Make an animation of the solution using Plot and ListAnimate.
In[51]:=
Click for copyable input
Out[52]=
A point is shown in red to demonstrate that the amplitude (pressure) of the time harmonic wave at a particular value of oscillates sinusoidally with the same frequency as the incoming wave.

Time Inharmonic Waves

Time inharmonic waves have arbitrary time variation, and may be represented by:
where is an arbitrary function, and its argument has the form in order to satisfy the wave equation (12). Again, the choice of defines the direction of travel of the wave. Setting makes the wave propagate in the negative direction and in the positive direction.
Explore a time inharmonic wave with adjustable speed of sound .
In[53]:=
Click for copyable input
Out[53]=PlayAnimation
There are three model parameters in the animation: a time slider, a direction switch and speed of sound selector. The wave propagation can be observed by playing the time slider.
Next, we build a 1D time inharmonic wave model using the wave equation (13) and solve it numerically using NDSolveValue.
The PDE for the sound pressure field is given by:
In[54]:=
Click for copyable input
Solve the PDE with NDSolveValue.
In[55]:=
Click for copyable input
Make an animation of the solution using Plot and ListAnimate.
In[56]:=
Click for copyable input
Out[58]=
In[59]:=
Click for copyable input
In the animation is shown that a point on the time inharmonic wave has arbitrary time variation.
Source Types
The wave equation model (14) contains two types of time harmonic sources: monopole and dipole sources. The monopole source models a point source that radiates sound equally in all directions. On the other hand, a dipole source is made up of two identical monopole sources that are in opposite phase and separated by a small distance compared to the wavelength of sound .

Monopole Sources

The monopole source from the wave equation (15) 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 [16]. The volume of sound medium is then displaced by the monopole source at a rate , which is known as the monopole source strength.
To make use of a monopole source, the monopole source is located at . Then the monopole source term may be written as
where is a Dirac delta function at the source location .
The Dirac delta function, however, poses a problem in numerical simulations as it can not be resolved in the discretized spatial domain. This is because the Dirac delta function is singular at the source location . Hence, an approximation to the Dirac delta function is needed. The process of approximating the Dirac delta function is called regularization.
There are various regularized delta functions available [17,18]. 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[60]:=
Click for copyable input
Inspect the regularized delta function centered at . Note that the width of support equals to .
In[61]:=
Click for copyable input
Out[61]=PlayAnimation
Modifying will change the width of the regularized delta function, however, for the spatial integral is always 1.
In[62]:=
Click for copyable input
Out[62]=
Next, we build a 1D time inharmonic wave model using the wave equation (19) and solve it numerically using NDSolveValue.
Assume the monopole source strength is sinusoidal with the frequency .
In[63]:=
Click for copyable input
Determine appropriate mesh spacing using the function introduced in Model Parameter Setup.
In[64]:=
Click for copyable input
Out[64]=
Define the monopole source as a function of , , and the location . Here we use the mesh spacing as the regularization parameter .
In[65]:=
Click for copyable input
Visualize the monopole source centered at .
In[66]:=
Click for copyable input
Out[69]=
Note that the amplitude in the animation depends on the selection of the regularized delta function , and should not be confused with the monopole source strength .
Set up a 1D acoustic model with the monopole source term .
In[70]:=
Click for copyable input
The PDE for the sound pressure field is given by:
In[71]:=
Click for copyable input
Solve the PDE with NDSolveValue.
In[72]:=
Click for copyable input
The analytical sound pressure amplitude of a monopole source is given by [20]:
Calculate the analytical sound pressure amplitude:
In[73]:=
Click for copyable input
Out[73]=
Visually verify the result.
In[74]:=
Click for copyable input
Out[76]=
It can be seen that the resulting wave has sinusoidal time dependence and that the monopole source is a point source at that radiates sound isotropically.

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 . Unlike a single monopole source, there is no net volume displacement by a dipole source. As one of the two monopoles expands the second contracts. A further difference between a monopole and a dipole is that a dipole source does not radiate sound isotropically.
To illustrate that a dipole is a combination of two monopole sources we model a dipole source at made from two monopoles separated by a distance . Based on the regularized delta function we choose in (21), the separating distance is automatically set as .
Set the separation distance of two monopoles.
In[77]:=
Click for copyable input
Out[77]=
In this case the regularization parameter needs to be small enough so that two monopoles do not overlap, since overlapping monopoles do not correctly model the dipole.
Inspect two monopoles with the regularization parameter ranging from to .
In[78]:=
Click for copyable input
Out[78]=PlayAnimation
Note that two monopoles begin to overlap when . The corresponding superposition (blue line) will no loner satisfy the integration rule and can not correctly model a dipole source.
Use as the regularization parameter to ensure that the two monopoles do not overlap. The negative sign of indicates that the two monopoles are in the opposite phase.
In[79]:=
Click for copyable input
The monopole source term in the wave equation (22) now becomes the summation of the two monopoles: .
Set up the monopole source term .
In[81]:=
Click for copyable input
Visualize the dipole source created from the two monopole sources.
In[82]:=
Click for copyable input
Out[85]=
The amplitude in the animation depends on the selection of the regularized delta function , and should not be confused with the monopole source strength .
Set up a 1D acoustic model using a pair of monopole sources.
In[86]:=
Click for copyable input
The PDE for the sound pressure field is given by:
In[87]:=
Click for copyable input
Solve the PDE with NDSolveValue.
In[88]:=
Click for copyable input
The analytical sound pressure amplitude of a dipole source is given by [23]:
Calculate the analytical sound pressure amplitude:
In[89]:=
Click for copyable input
Out[89]=
Visually verify the result.
In[90]:=
Click for copyable input
Out[92]=
Compute the difference of the amplitudes of a pair of monopole sources at and the analytical value.
In[93]:=
Click for copyable input
Out[93]=
The dipole source does not radiate sound isotropically. The resultant sound waves are sinusoidal but in opposite phase.
In the wave equation , the dipole source is equivalent to the two monopole sources just demonstrated. The dipole source in the wave equation (24) is written as (25), where represents the dipole source strength.
The dipole source strength is equal to the product of the monopole pair strength , the medium density and the separation distance .
Just as for monopoles the Dirac delta function needs to be regularized. Note that the choice for that we previously defined has a continuous first derivative so it can also be used to model a dipole source that appears in the divergence term in the wave equation (26).
Set up the 1D dipole source .
In[94]:=
Click for copyable input
The amplitude in the animation depends on the selection of the regularized delta function , and should not be confused with the dipole source strength
Visualize the dipole source .
In[96]:=
Click for copyable input
Out[99]=
The resulting sound field is not isotropic, but is symmetric about the axis connecting the two monopoles.
Set up the acoustic model with a dipole source term .
In[100]:=
Click for copyable input
The PDE for the sound pressure field is given by:
In[101]:=
Click for copyable input
Solve the PDE with NDSolveValue.
In[102]:=
Click for copyable input
Analytically the dipole source is identical to the pair monopoles with a corresponding source strength. However, since we applied the regularized delta function to approximate the true delta function in a discretized domain, there is a difference between the single dipole model and the model made up from a pair of monopoles. This difference will only be present within the support of the regularized delta function .
Show the difference between the dipole source and the pair monopoles at .
In[292]:=
Click for copyable input
Out[292]=
Show the difference between the dipole model and pair monopoles model (blue) superimposed on the wave scaled to fit on the same vertical scale (light gray).
In[104]:=
Click for copyable input
Out[106]=
As expected, the dipole model is consistent to the pair monopoles model outside the support of the regularized delta function (), which means a dipole is a combination of a pair of monopole sources , and can be directly modeled by a dipole source term .
The Wave Equation as a System of First Order Equations
The following section introduces an alternative way to model sound propagation, which makes use of a system of first order differential equations. Sometimes using a system of first order equations is the only way to solve an acoustic model, for example when perfectly matched layers (PML) are to be used. Making use of a first order reformulation of the the PDE might also increase the time and memory efficiency during the solution process.
Assume that sound is transmitted through a medium by longitudinal waves. This means that particles move back and forth in the propagation direction. The movement of these particles is denoted by the sound particle displacement , and its time derivative gives the sound particle velocity .
Note that the sound particle velocity is not the same as the speed of sound , that is the speed at which longitudinal waves propagate.
The equations of mass and momentum conservation can be formulated by Newton's law and the ideal gas law [27]:
The coupled first order wave equations are the combination of (28) and (29).
As an example consider a 1D model of time harmonic waves.
Set up a 1D source free acoustic model using the coupled first order wave equations.
In[107]:=
Click for copyable input
Note that when we build the model with the coupled first order wave equations, extra initial and boundary conditions are required for the sound particle velocity
To set up the extra initial and boundary condition, the specific impedance should be defined first. The specific impedance denotes the ratio of the sound pressure to the sound particle velocity as the wave moves through the medium, which is written as:
For a plane wave the specific impedance can be shown to be equal to the product of the mediums density and the speed of sound [30]. A detailed explanation of impedance can be found in the section entitled Impedance Boundary Conditions.
The extra boundary condition may be formulated with the specific impedance using (31)
Set the impedance of air for 1D plane wave.
In[108]:=
Click for copyable input
In this case the initial conditions are set to zero. Boundary conditions are enforced so that a sinusoidal pressure source is placed at the left end of the domain.
Set the initial conditions and boundary conditions.
In[109]:=
Click for copyable input
The PDE for the sound pressure field is given by:
In[111]:=
Click for copyable input
Solve the system of PDEs and monitor time/memory usage.
In[112]:=
Click for copyable input
Make an animation of the solution using Plot and ListAnimate.
In[114]:=
Click for copyable input
Out[116]=
It is instructive to compare the result with the one computed with a second order wave equation (32).
The same sinusoidal pressure source is placed at the left boundary.
In[117]:=
Click for copyable input
The PDE for the sound pressure field is given by:
In[118]:=
Click for copyable input
Solve the second order PDE and monitor the time/memory usage.
In[119]:=
Click for copyable input
Set up the analytical solution.
In[121]:=
Click for copyable input
Compare the error of the coupled first order PDEs and second order PDE at .
In[122]:=
Click for copyable input
Out[122]=
Note that the error, memory and time consumption are greatly reduced when we use the first order system of PDEs. This is very helpful when we are solving a large scale acoustic models.
Sound Propagation in Lossy Media
To model sound propagation in lossy media such as porous materials, the wave equation (33) needs to be modified in order to accommodate for an attenuation of the sound waves [34]. The effective bulk modulus and an additional first order time derivative on are added on the left hand side.
Here the attenuation term and the effective bulk modulus account for the effects of porosity, flow resistivity, together with any deviation from adiabatic compressibility. The porosity is defined as the ratio of the void spaces to the total volume of a medium, and the flow resistivity is defined as the ratio of the pressure gradient to the sound particle velocity . For common sound absorbent materials values typically lie in the range of to .
In the lossy media may deviate from the adiabatic bulk modulus due to the dissipation of the sound. The value of varies with frequency and the flow resistivity , and can be approximated by an empirical formula from Delany and Bazley [35]. Note that equation (36) simplifies to the lossless wave equation (37) when , , and .
Since the effective bulk modulus depends on a specific frequency , the lossy media model (38) is constrained to simulate single frequency waves only, that is, a time-harmonic wave. To model a time-inharmonic wave that composes of multiple frequencies, the wave should be decomposed by a Fourier transform before analysis.
As an example we consider a harmonic sound wave at frequency , which propagates 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[123]:=
Click for copyable input
Calculate the effective bulk modulus .
In[124]:=
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 to accommodate for the modified wave equation (39) [40].
Set up an absorbing boundary condition for the lossy media model.
In[125]:=
Click for copyable input
A harmonic sound source at frequency is added on the left end using DirichletCondition.
In[126]:=
Click for copyable input
The PDE for the sound pressure field is given by:
In[128]:=
Click for copyable input
Solve the PDE with NDSolveValue.
In[129]:=
Click for copyable input
For comparison set up an analytical time harmonic wave propagating in a lossless media.
In[130]:=
Click for copyable input
Compare the wave propagation in a lossy media (blue) and in a lossless media (gray).
In[131]:=
Click for copyable input
Out[133]=
In the lossy media the sound wave propagates in a lower speed while being attenuated within the domain.
Acoustic Boundary Conditions
The most common boundary conditions in acoustics can be modeled with DirichletCondition and NeumannValue, and can be categorized in the following three types:
Under these three types, the following boundary conditions are introduced:
The following section describes several physical boundaries common in acoustics and how they can be modeled with the use of DirichletCondition and NeumannValue. For this purpose the boundary condition currently discussed is always on the right hand side of the simulation domain. In some examples additional boundary conditions are on the left hand side to better demonstrate the behaviour of the boundary condition on the right hand side.
Impedance Boundary Conditions
When a sound wave propagates from one medium into another, part of the sound wave will be reflected at the interface and part will be transmitted across the interface. To illustrate this behaviour a tube filled with two sound mediums is considered in the next example.
For a tube aligned with the axis and a diameter much smaller than the sound wave length , it is then reasonable to assume that all points in the , plane have the same sound pressure value and the same sound particle velocity . This wave type is called plane wave and can be treated as a one-dimensional wave.
To better illustrate an impedance boundary condition a Gaussian pulses with standard deviation is used in this example.
Set the initial conditions of the sound pressure field.
In[134]:=
Click for copyable input
The region is split at into two media. The left medium is cold air and the right medium is hot air .
Specify speeds of sound and densities with ThermodynamicData.
In[137]:=
Click for copyable input
In this example the simulation domain is enlarged a little. This way the wave will not reach the right boundary before the end of the simulation time.
Define an extended domain from to that will include two media.
In[141]:=
Click for copyable input
The PDE for the sound pressure field is given by:
In[142]:=
Click for copyable input
Solve the PDE with NDSolveValue.
In[143]:=
Click for copyable input
Visualize the solution and inspect the region near the interface at .
In[144]:=
Click for copyable input
Out[146]=
The sound pressure field animation shows a sound wave moving from left to right, and two mediums form an interface in the middle at . At the interface the wave is reflected with the same speed but opposite polarity. The remaining part of the wave transmits across the interface with higher speed but smaller amplitude. The wave amplitude change depends on the material properties of both regions. The reflected wave behaves a bit like a sound soft boundary. This is because of the difference in density and . Changing the densities such that will result in a different reflection. This can be checked by simply swapping , with , .
One property that is used to formulate the relation between the incident and the reflected and the 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:
For a plane wave the specific impedance can be shown to be equal to the product of mediums density and the speed of sound [41].
For a sound wave propagating from one medium to another, the amplitude relations can be expressed in specific impedance as
where , and are the amplitude of the reflected, the incident and the transmitted wave, respectively. and are the specific impedance of the medium 1 and 2, respectively [42].
Compute the specific impedances and of this example.
In[147]:=
Click for copyable input
Out[147]=
Compute reflected and transmitted amplitudes for an incident amplitude of .
In[148]:=
Click for copyable input
Out[148]=
With these parameters the solution can be verified.
Show the solution at with the predicted transmitted and reflected amplitudes.
In[149]:=
Click for copyable input
Out[149]=
Compute the error in the amplitude of the reflected wave at .
In[150]:=
Click for copyable input
Out[150]=
Compute the error in the amplitude of the transmitted wave at .
In[151]:=
Click for copyable input
Out[151]=
If one is not interested in the solution in, say, medium then there is no need to include that region in the solution. The medium interface can be modeled by specifying the impedance of medium at the boundary at . This is called an impedance boundary condition. An impedance boundary condition can be realized with a generalized Neumann boundary condition which is also known as a Robin boundary condition.
Substituting the definition of impedance (43) into the momentum conservation equation (44) , then the impedance boundary condition can be expressed in the boundary impedance as
Rebuild the example with known boundary impedance , so the NeumannValue is set to . Note that the coefficient in needs to match the NeumannValue specification . Here and are 0 and matches and this is then consistent with the wave equation (45) we use. Also, the medium boundary is at the right end where . The meaning of time derivatives in NeumannValue is explained in more detail in the section entitled Neumann Value on Time Derivative.
Set up the impedance boundary condition with NeumannValue.
In[152]:=
Click for copyable input
The PDE for the sound pressure field is given by:
In[153]:=
Click for copyable input
Solve the sound pressure field only in the region of medium , not including the part with medium 2 included in the used above.
In[154]:=
Click for copyable input
Make an animation of the solution using Plot and ListAnimate.
In[155]:=
Click for copyable input
Out[157]=
The animation shows the same result as the previous example just with a smaller simulation domain.
Visually verify the solution at .
In[158]:=
Click for copyable input
Out[158]=
Compute the error in the amplitude of the reflected wave at .
In[159]:=
Click for copyable input
Out[159]=
Using an impedance boundary condition in principal allows a simulation process to become more efficient since there is no longer the need to solve the equation in both media.
Conceptually, the impedance boundary condition
is a generalization of three types of boundary conditions.
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 boundaries at infinity. An absorbing boundary condition (ABC) works by absorbing an incoming wave and thus make the model behave as if it had infinite extent. ABCs are not the only way to model simulation domains with infinite extent. As an alternative approach to ABCs perfectly matched layers (PML) may be used.
Since an absorbing boundary is a continuation of the computational domain, the boundary impedance is equal to the specific impedance of the domain. That is, an ABC is a special case of an impedance boundary condition that make use of the time derivative of the dependent variable within the NeumannValue. The meaning of time derivatives in NeumannValue is explained in more detail in the section entitled Neumann Value on Time Derivative.
Recall that for a plane wave the specific impedance is . Substituting into the impedance boundary condition (47), then the absorbing boundary condition can be written as
To accommodate for different kinds of waves, equation (48) can be generalized with an additional term : [49]
For a plane wave , for a cylindrical wave and for a spherical wave are used.
Here is the distance from the wave origin to the absorbing boundary.
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 an absorbing boundary condition at the right end.
Inspect the setting of a plane wave absorbing boundary condition on the right end.
In[160]:=
Click for copyable input
Out[160]=
The PDE for the sound pressure field is given by:
In[161]:=
Click for copyable input
Solve the PDE with NDSolveValue.
In[162]:=
Click for copyable input
Make an animation of the solution using Plot and ListAnimate.
In[163]:=
Click for copyable input
Out[165]=
Note that the incoming wave is absorbed at the right hand boundary as if the simulation domain had infinite extent.
One quality characteristic of a good absorbing boundary condition is the amount of how much an incoming wave is reflected back into the simulation domain. Ideally nothing is reflected back.
Set up the analytical solution.
In[166]:=
Click for copyable input
Inspect the error (blue) propagation at various times of the scaled wave (light gray).
In[168]:=
Click for copyable input
Out[170]=
The numerical error is due to the discretization of the finite element method. Note that the error didnt increase as the wave reached the right end, which means there is very little to no reflection at the absorbing boundary in this 1D case.
However, the absorbing boundary condition can only absorb waves exactly at normal incidence. In 2D or 3D cases where waves can strike a boundary at any possible angle, a numerical reflection will occur at the absorbing boundary. This is illustrated in a simulation of a cylindrical wave propagation in a 2D domain with absorbing boundaries.
Set a 2 dimensional, time dependent wave equation without sources in inactive form.
In[171]:=
Click for copyable input
Define the 2D domain.
In[172]:=
Click for copyable input
Set up an absorbing boundary condition for a cylindrical wave at all four boundaries.
In[174]:=
Click for copyable input
Set the initial conditions of the sound pressure field. A 2D Gaussian pulse with the standard deviation is used.
In[176]:=
Click for copyable input
The PDE for the sound pressure field is given by:
In[179]:=
Click for copyable input
Solve the PDE with NDSolveValue and monitor time/memory usage.
In[180]:=
Click for copyable input
To visualize the wave reflection, we set up a legend bar and ContourPlot options scaling from to .
In[182]:=
Click for copyable input
Visualize the truncated solution with the vertical PlotRange set sufficiently small to see the reflections.
In[184]:=
Click for copyable input
Out[187]=
A contour plot shows the size of the reflected wave at .
In[188]:=
Click for copyable input
Out[188]=
Note that the reflected wave at the four corners where the incident angles are not normal to boundaries. With a different geometry the reflections may be somewhat smaller.
To reduce the wave reflection in 2D or 3D, a perfectly matched layer (PML) can be applied as an alternative method to model an absorbing boundary. A PML is designed such that waves that impact boundaries in a non-normal way are absorbed within it. Further details on PML can be found in the section entitled Perfectly Matched Layer (PML).
Sound Hard Boundary Conditions - Walls
The sound particle velocity on the boundary corresponds to the speed of a mechanical vibration. For sound hard boundaries, the sound particle displacement at the boundary is fixed to zero . As a consequence the sound particle velocity component normal to the boundary is also equal to zero.
At a sound hard boundary the impedance approaches infinity since the sound velocity . That is, the sound hard boundary condition is a special case of an impedance boundary condition:
Derived from Newton's law [50], the momentum conservation equation (51) relates the sound pressure with the sound particle velocity .
Combining (52) and (53), it is found that a sound hard boundary condition also implies that the partial derivative of the sound pressure in direction normal to the boundary is zero.
As an example we look at a tube with one end closed.
To model a sound hard boundary condition the NeumannValue is set to at the closed end of the tube at .
In[189]:=
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 sound hard boundary are the default boundary condition used if no boundary condition is specified at a given boundary.
In[190]:=
Click for copyable input
The PDE for the sound pressure field is given by:
In[191]:=
Click for copyable input
Solve the PDE with NDSolveValue.
In[192]:=
Click for copyable input
Make an animation of the solution using Plot and ListAnimate.
In[193]:=
Click for copyable input
Out[195]=
This sound pressure field animation shows a sound wave moving from left to right towards a sound hard boundary. Conceptually a hard sound boundary is equivalent to a closed end of a tube. Since no forward motion is possible, the pressure accumulates and reaches its maximum at the double of the incident pressure. The pressure wave then starts moving backward. Note that the wave was reflected with the same speed, amplitude and polarity as the incident wave.
Normal Velocity Boundary Conditions
A boundary where the sound particle velocity normal to the boundary is specified and not equal to zero is called a normal velocity boundary. A familiar example of a normal velocity boundary is a vibrating mechanical structure that produces sound waves into the surrounding medium.
The negative sign is added in front of to indicate that the sound particle velocity is specified opposite to the outward normal . By doing so we ensure that the sound pressure will always be in phase with the particle velocity on a boundary.
Recall that the momentum conservation equation (54) relates the sound pressure with the sound particle velocity .
Inserting (55) into (56), then the normal velocity boundary condition can be written as
In the following example a vibration is introduced at the right hand boundary and vibrates with a known velocity . The sound field can be calculated by using NeumannValue as shown below.
Specify the vibration velocity.
In[196]:=
Click for copyable input
Set up the normal velocity condition at the right end with NeumannValue.
In[197]:=
Click for copyable input
The PDE for the sound pressure field is given by:
In[198]:=
Click for copyable input
Solve the PDE with NDSolveValue.
In[199]:=
Click for copyable input
Make an animation of the solution using Plot and ListAnimate.
In[200]:=
Click for copyable input
Out[202]=
The simulation begins with initial conditions are rest and the normal velocity boundary condition generates a sinusoidal wave on the right hand boundary that propagates to the left. The small wiggles ahead of the wave are artifacts due to the spatial discretization. This effect can be reduced by setting a higher mesh resolution in AcousticOptions.
Note that this animation is quite similar with the example shown in Pressure Source Boundary Conditions. That is because the sound pressure and particle velocity have linear dependence with each other given by , where is the specific impedance, so the perturbation on the sound pressure will also cause the perturbation on the sound particle velocity. When the specific impedance is known a normal velocity boundary can be replaced by a pressure source boundary.
Sound Soft Boundary Conditions
On a sound soft boundary the medium pressure at the boundary is set equal to an ambient reference pressure such that . This type of boundary is also called sound free boundary since the sound particles are free to move at the boundary, which means that the sound particle velocity is not required to be zero. The momentum conservation equation relates the sound particle displacement and sound pressure . While , is not necessarily equal to zero.
At a sound soft boundary the impedance since the sound pressure . That is, the sound soft boundary condition is a special case of impedance boundary conditions when the boundary impedance is set to zero:
To illustrate this behaviour a tube with one end open is considered in the next example. The open-ended side is treated as a sound soft boundary since there is no constraint to limit the sound particles movement.
To model the sound soft boundary condition a DirichletCondition is set to at the right open end of the tube.
In[203]:=
Click for copyable input
The PDE for the sound pressure field is given by:
In[204]:=
Click for copyable input
Solve the PDE with NDSolveValue.
In[205]:=
Click for copyable input
Make an animation of the solution using Plot and ListAnimate.
In[206]:=
Click for copyable input
Out[208]=
The sound pressure field animation shows a sound wave moving from left to right. At the right end a sound soft boundary is positioned. Conceptually a sound soft boundary is equivalent to an open end of a tube. As the wave reaches the open end its pressure quickly drops to the ambient pressure, which creates a suction behind it and forms a negative pressure pulse, and then the negative pressure pulse starts moving backward as the reflected wave. Noted that the wave is reflected with the same speed and amplitude as the incident wave, but the phase angle is changed by .
Pressure Source Boundary Conditions
We speak of a pressure source boundary condition when a nonzero, possibly time dependent sound pressure is specified at a boundary. Both DirichletCondition and NeumannValue can be used to specify a pressure source boundary condition.

Dirichlet Model

To model, for example, a vibrating panel that produces a time dependent sinusoidal pressure wave, a pressure source can be set up as a time dependent sinusoidal pressure at the right end. And an absorbing boundary condition is added at the left end to avoid wave reflection.
Set up the pressure source boundary with DirichletCondition.
In[209]:=
Click for copyable input
The PDE for the sound pressure field is given by:
In[211]:=
Click for copyable input
Solve the PDE with NDSolveValue.
In[212]:=
Click for copyable input
Make an animation of the solution using Plot and ListAnimate.
In[213]:=
Click for copyable input
Out[215]=
The simulation begins with an undisturbed domain and the transient pressure variation generates a sinusoidal wave on the right hand boundary that propagates to the left.
Note that a pressure source generates waves normal to the boundary while a monopole source radiates waves isotropically. So for 1D cases, a pressure sources on the boundary can be modeled by replacing the DirichletCondition with a equivalent monopole source on the right hand side of the wave equation (57). Monopole sources are explained in the section entitled Monopole Sources.

Neumann Model

A pressure source can also be modeled with a NeumannValue. Doing so has the advantage that on the same portion of the boundary the pressure source can be combined with other Neumann values based boundary condition models. In other words hybrid boundary conditions can be created by combining several Neumann values on the same part of the boundary. In this way, for example, a pressure source can be combined with an absorbing boundary to form a radiation boundary condition.
The NeumannValue setting for a pressure source can be derived from the normal velocity boundary conditions:
Recall that for a plane wave the specific impedance . Substituting into (58) the equation can be expressed in terms of the specified sound pressure as
As an example we model a pressure source with the NeumannValue setting (59).
Set up the pressure source boundary with NeumannValue.
In[216]:=
Click for copyable input
The PDE for the sound pressure field is given by:
In[217]:=
Click for copyable input
Solve the PDE with NDSolveValue.
In[218]:=
Click for copyable input
Make an animation of the solution using Plot and ListAnimate.
In[219]:=
Click for copyable input
Out[221]=
Set up the analytical solution.
In[222]:=
Click for copyable input
Compare the difference between the Dirichlet model and the Neumann model at .
In[223]:=
Click for copyable input
Out[223]=
Note that the Dirichlet mode is more accurate because Dirichlet conditions enforce a specific values at the boundary nodes. This means that modeling a pressure source with a Neumann value should only be done if there is a need to combine that pressure source with a other Neumann value based boundary conditions.
Radiation Boundary Conditions
Conceptually a radiation boundary is a hybrid boundary condition that combines the properties of a pressure source boundary and an absorbing boundary. A nonzero, possibly time dependent 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.
To derive the radiation boundary condition a sound pressure field is broken down into the incoming wave and outgoing wave as
One way to produce an incoming wave is by using the pressure source boundary condition with a specified sound pressure at the boundary.
An absorbing boundary condition is applied to ensure that the outgoing wave leaves the boundary without constraint.
Adding up equations (60) and (61) and inserting (62) gives
Recall that on the boundary the amplitude of incoming wave is specified by . Substituting this relation into the right hand side of (63) yields the radiation boundary condition
To illustrate this behaviour a semi-infinite tube with one end closed is considered in the next example. A harmonic sound wave enters the computational domain from the right end. The wave then travels through the tube. The left end is then treated as a sound hard boundary. When the wave comes back to the right end it is absorbed. This process of being a pressure source and absorbing the reflected wave is the radiation boundary condition.
Set up the radiation boundary with NeumannValue at the right end. By default a sound hard boundary condition is implicitly used at the left end.
In[224]:=
Click for copyable input
The PDE for the sound pressure field is given by:
In[226]:=
Click for copyable input
Solve the PDE with NDSolveValue. To show the effect of the radiation boundary the simulation time is extended a little.
In[227]:=
Click for copyable input
Make an animation of the solution using Plot and ListAnimate.
In[229]:=
Click for copyable input
Out[231]=
The animation shows a left-traveling wave that is produced by the radiation boundary. The wave is then reflected at the left end and superimposes with the incident wave. As the reflected wave starts leaving through the radiation boundary the sound pressure field reaches its dynamic steady state. The resulting wave neither moves right nor left but simply oscillates in time.
This type of wave is called standing wave, and will be explained in more detail in the separate tutorial entitled Acoustics in the Frequency Domain.
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. A PML, however, is not a boundary condition. The PML is an artificial additional absorbing layer that enlarges the original simulation grid. As a wave enters the absorbing layer, it is attenuated by an absorption and decays exponentially [64]. The key property of a PML that distinguishes it from an ordinary absorbing boundary condition is that it is designed so that waves propagate from a non-PML region to a PML and do not reflect at the interface in any incident angle. That is, a PML acts like a non-reflecting absorbing material. Further details on PML can be found in the appendix section entitled Perfect Matched Layer Derivation.
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.
It is better to split the wave equation into a coupled differential equations with first-order only spatial derivatives. Then the entire process of generating a PML can be performed by a single transformation (65) of the original wave equation.
The coupled first-order wave equation after the PML coordinate transformation is
The absorption coefficient is a 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.
A PML interface is only reflectionless when we are solving the exact, continuous wave equation. NDSolve, however, is a numerical method that discretizes the simulation domain and solves an approximation to the wave equation instead. In this numerically approximate case the wave will show reflections at the PML interface. To reduce the amount of reflection a finer grid or a smaller attenuation rate (smaller ) should be used.
For an example we consider a simulation of a 1D domain from to using a computational domain that from to with a PML at the right from to .
Define the 1D domain and the PML width.
In[232]:=
Click for copyable input
Set up the coupled first-order wave equation with PML.
In[234]:=
Click for copyable input
Take air as the sound medium, and for an one dimensional plane wave the specific impedance is [66].
Compute the specific impedance .
In[235]:=
Click for copyable input
Out[235]=
The initial conditions are set to zero. Boundary conditions are enforced so that a sinusoidal pressure source is placed at the left end of the domain. Recall that the sound particle velocity is related to the sound pressure by . Details of the setup as coupled first-order wave equations can be found in the section entitled The Wave Equation as a System of First Order Equations.
Set the initial conditions and boundary conditions.
In[236]:=
Click for copyable input
Specify and set the absorption coefficient to increase linearly within the PML layer.
In[238]:=
Click for copyable input
Inspect the constructed absorption coefficient .
In[240]:=
Click for copyable input
Out[240]=
As seen in the original non-PML region.
The PDE for the sound pressure field is given by:
In[241]:=
Click for copyable input
Solve the PDE with NDSolveValue. To show the effect of PML the simulation time is extended a little.
In[242]:=
Click for copyable input
Visualize the result including the PML region.
In[244]:=
Click for copyable input
Out[247]=
A thinner PML requires a greater value of to get the same attenuation. Larger values of generally give larger numerical reflection. To inspect the tradeoff between PML width and the attenuation rate, this PML model is solved with different parameters as follows
Experiment with different values for the PML width and the attenuation coefficient .
In[248]:=
Click for copyable input
Create a helper function fun that takes the PML width and as arguments.
In[250]:=
Click for copyable input
Solve the PDE with different PML parameters.
In[251]:=
Click for copyable input
Show the results of different PML settings.
In[252]:=
Click for copyable input
Out[254]=
For small a wider PML is required to attenuate the wave, however, for large the numerical reflection is much more apparent.
The main advantage of PML is to absorb incident waves isotropically. To show this consider a 2D example of a unit square with PML at four boundaries. The simulation of a 2D domain is using a computational domain with the PML at and .
Define a 2D domain with PML at four boundaries.
In[255]:=
Click for copyable input
Specify and set the absorption coefficients and that act on each dimensions.
In[259]:=
Click for copyable input
To implement the PML on a 2D wave equation two auxiliary field variables and are introduced [67]. The equation becomes:
Here and are sound particle velocity components in space.
Set up a 2 dimensional, first-order wave equation with PML.
In[262]:=
Click for copyable input
Set the initial conditions of the sound pressure field. A 2D Gaussian pulse with the standard deviation is used.
In[263]:=
Click for copyable input
The PDE for the sound pressure field is given by:
In[266]:=
Click for copyable input
Solve the PDE with NDSolveValue and monitor time/memory usage.
In[267]:=
Click for copyable input
To visualize the wave reflection, we set up the legend bar and the ContourPlot options scaling from to .
In[269]:=
Click for copyable input
Visualize the truncated solution in the non-PML region with the PlotRange set to .
In[271]:=
Click for copyable input
Out[274]=
Compared to the previous 2D case shown in the section Absorbing Boundary Conditions, the wave is now better absorbed at the corners with PML at the cost of extra time/memory usage.
Appendix
Perfectly Matched Layer Derivation
The perfectly matched layer (PML) is an artificial absorbing layer that is placed adjacent to the edges of the grid. As a wave enters the absorbing layer, it is attenuated by the absorption and decays exponentially [68]. The key property of a PML that distinguishes it from ordinary absorbing material is that it is designed so that waves propagate from a non-PML medium to PML do not reflect at the interface. That is, a PML acts like a non-reflected absorbing material.

Attenuation using an Artificial Complex Dimension

Recall that a time harmonic wave traveling in the positive x-direction is expressed as
The concept of a PML is to perform a coordinate transformation that all the terms in wave expression (69) are replaced by a complex variable .
where the original serves as the real part of , and is an artificial imaginary part in an absorbing layer.
This way a wave will decay exponentially in PML but remain unchanged otherwise. The PML attenuation rate is controlled by the slope of the artificial imaginary part . For a thinner PML a greater attenuation rate is required, that is, a larger value has to be assigned. However, using a very large value can increase the numerical reflection of the PML, so there is a tradeoff between the PML width and its associated quality and the simulation speed.
For an example we consider a simulation of a 1D domain from to using a computational domain that from to with a PML at the right from to . To achieve this we add a linearly growing imaginary part with the slope in the PML.
Define the 1D domain and the PML width.
In[275]:=
Click for copyable input
Manually chosen slope of the artificial imaginary part .
In[277]:=
Click for copyable input
Show in the complex plane.
In[279]:=
Click for copyable input
Out[279]=
With the PML coordinate transformation (70) the wave expression (71) transforms to
Since there is an extra term in the new wave expression (72), the wave with positive wave number will exponentially decay in the absorbing layer where .
Calculate the wave number in air.
In[280]:=
Click for copyable input
Inspect the sound pressure at when PML is placed on .
In[281]:=
Click for copyable input
Out[281]=
Note that the wave is unchanged outside the PML () where .
To increase the attenuation rate in PML, we can pick a greater slope value for the artificial imaginary part .
In[282]:=
Click for copyable input
In[284]:=
Click for copyable input
Out[284]=
Now the wave is attenuated rapidly and decays to nearly zero before . In other words, the greater attenuation rate (larger value) allows us to use a thinner PML.
Inspect the wave with time propagation.
In[285]:=
Click for copyable input
Out[287]=
Note that there is no reflection at the PML interface (), so it acts like a non-reflecting absorbing material.

PML Coordinate Transformations on the Coupled 1st Order Wave Equation

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. Recall that for the coordinate transformation all terms in the wave equation (73) are replaced by the complex valued variable defined as:
In order to make the PML attenuation rate independent of the wave frequency [74], we set the slope to be
Here is the absorption coefficient defined by users. To minimize the numerical reflection in the discretized domain, is set to increase linearly within the PML layer from to . The coordinate transformation becomes
or
However, the original wave equation (75) contains the second order spatial derivative. So it is more convenient to split the wave equation into coupled differential equations with only first-order spatial derivative, then the entire process of PML can be performed by a single transformation (76) of the original wave equations.
Recall that the coupled first-order wave equation is shown as
So in the 1D case the coupled wave equation becomes
Now we perform the PML coordinate transformation (77) on (78), and apply the harmonic wave relation and . The equations become
Multiplying both sides by yields
Apply the harmonic wave relation again to turn the equations back into time domain form. Finally, the coupled first-order wave equations with PML become
NeumannValues on Time Derivatives
An absorbing boundary condition can be modeled with a NeumannValue on the time derivative. This is best understood when looking at a simplified wave equation.
(79) can be written as a system of equations; and in fact this transformation is done internally in NDSolve. For the transformation we introduce the following auxiliary equation:
If we now substitute (80) into (81) and use (82) we obtain an equivalent system of first order equations for the original wave equation (83).
Set up initial conditions for a wave equation.
In[288]:=
Click for copyable input
Set up the equation.
In[289]:=
Click for copyable input
Next, we perform the substitution from (84) and modify the initial conditions.
Set up the modified initial conditions.
In[290]:=
Click for copyable input
Set up the system of first order equations.
In[291]:=
Click for copyable input
The important point is to note how the time derivative in the Neumann value changed by applying (85) to the auxiliary dependent variable. In this case the time derivative in the Neumann value is replaced by to become a Neumann value on the auxiliary variable.
Nomenclature
The following table summarizes symbols names used and their meaning in this tutorial.
    

SymbolDescriptionUnit
ρdensity of a medium[kg/m3]
cspeed of sound in a medium[m/s]
psound pressure[Pa]
pmaxmaximum value of sound pressure[Pa]
sound pressure amplitude function[Pa]
specified boundary pressure [Pa]
ttime[s]
tendsimulation end time[s]
Xposition vector[m]
sdirection switchN/A
Foptional dipole source[N/m3]
dipole source strength[N/m3]
Qoptional monopole source[1/s2]
monopole source strength[1/s2]
Nullseparation 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]
αattenuation factor[m2/(s·N)]
ϕporosityN/A
VVvoid volume[m3]
VTtotal volume[m3]
Rfflow resistivity[kg/(m3·s)]
βstandard deviation of a Gaussian pulse[m]
ζsound particle displacement[m]
Nullsound particle velocity[m/s]
Nullspecified boundary velocity [m/s]
Ttemperature[K]
Zspecific impedance[Null]
Zbboundary impedance[Null]
Aramplitude of reflected wave[Pa]
Aiamplitude of incident wave[Pa]
Atamplitude of transmitted wave[Pa]
complex spatial variable of PML[m]
fpPML model parameter[m]
mfslope of PML artificial imaginary partN/A
σabsorbtion coefficient of PML[rad/(s·m)]
σmaxmaximum value of absorbtion coefficient[rad/(s·m)]
Nullarbitrary functionN/A

References

1.  Nolte, Bodo and Marburg, Steffen. Computational acoustics of noise propagation in fluids : finite and boundary element methods. Springer - Verlag, 2008.

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

3.  Russell, Daniel. Acoustics and Vibration Animations: Reflections from Impedance and the Standing Wave Ratio, The Pennsylvania State University, 2013, www.acs.psu.edu/drussell/Demos/SWR/SWR.html.

4.  Saksela, Kai. Finite element analysis of 2D acoustics. Kais thoughts, 2013, blog.kaistale.com/?p=574.

5.  Fahy, Frank. Foundations of Engineering Acoustics. Academic Press, 2001.

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

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

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

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

10.  Zeleny, Enrique. Acoustic Multipoles, Wolfram Demonstrations Project, 2013, http://demonstrations.wolfram.com/AcousticMultipoles/.

11.  Zeleny, Enrique. Longitudinal and Transverse Waves, Wolfram Demonstrations Project, http://demonstrations.wolfram.com/LongitudinalAndTransverseWaves/.

12.  Moore, Guy. Lecture Note: The Physics of Music - Reflection and Impedance, McGill University, 2006.

13.  MATTER Project. Full Width Half Maximum (FWHM), The University of Liverpool, Department of Physics, 2000.

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

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

16.  Maderuelo-Sanz et al. Acoustical performance of porous absorber made from recycled rubber and polyurethane resin. Latin American Journal of Solids and Structures, Vol.10, No.3, 2013.

17.  Turkel, E. and Yefet, A. Absorbing PML boundary layers for wave-like equations. Applied Numerical Mathematics, 27(4) 533-557, 1998.

18.  T. Cox, P. D'Antonio. Acoustic Absorbers and Diffusers: Theory, Design and Application. Spon Press, 2004.

19.  B. Engquist and A. Majda. Absorbing Boundary Conditions for the Numerical Simulation of Waves. Mathematics of Computation, Vol. 31, No. 139 629-651, 1977.