Acoustics in the Time Domain
| Contents | Perfectly Matched Layer |
| Introduction | Appendix |
| Wave Equation | Nomenclature |
| Acoustic Boundary Conditions | References |
Contents
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, 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".
Extended examples of sound system modeling can be found in the Model Collection.
Many of the animations of the simulation results shown in this notebook are generated with a call to Rasterize. This is to reduce the disk space that this notebook requires. The downside is that the visual quality of the animations will not be as crisp as without it. To obtain high-quality graphics, remove or comment out the call to Rasterize.
The symbols and corresponding units used throughout this tutorial are summarized in the Nomenclature section.
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.
Various sections in the documentation explain the use of inactive PDE operators. Please refer to "Numerical Solution of Partial Differential Equations".
The actual numerical material parameters can be put into the model using replacement:
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. A function is created to set the max cell measure for both time-harmonic and time-inharmonic waves.
Throughout this tutorial, a particular example frequency will be used for time-harmonic waves and a Gaussian pulse as an example of time-inharmonic waves.
The default resolution of the max edge length is set to 12 nodes per wavelength 
, which means that there will be at least 12 elements per wavelength 
in each direction of the wave propagation. 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.
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].
The following model parameters are used for the examples in this tutorial. These parameters define the simulation domain 
and the simulation end time 
.
If an acoustic wave begins with an undisturbed domain, the initial sound pressure 
and its derivative are set to zero.
In the examples of time-inharmonic waves, typically Gaussian pulses will be used, making use of a standard deviation of 
.
In many examples, an absorbing boundary condition (ABC) will be used to avoid wave reflection.
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 behavior, 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 
direction.
Combining (8) and (9), a traveling time-harmonic wave is expressed as:
and wave number 
. Assume that the amplitude 
: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 increases when selecting a higher-frequency 
. The number of phase changes per unit distance increases as the wave number 
increases. The wave propagation can be observed by playing with 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, a 1D time-harmonic wave model is built using the wave equation (11) and solved numerically using NDSolveValue.
will produce a time-harmonic wave. The source can be added using DirichletCondition.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.
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 with the time 
slider.
Next, a 1D time-inharmonic wave model is built using the wave equation (13) and solved numerically using NDSolveValue.
The animation shows 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, 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 
.
Modifying 
will change the width of the regularized delta function; however, for 
, the spatial integral is always 1:
Next, a 1D time-inharmonic wave model is built using the wave equation (19) and solved numerically using NDSolveValue.
. Here the regularization parameter is defined to be half of the mesh spacing: 
:
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 
.
The analytical sound pressure amplitude of a monopole source is given by [20]:
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, a dipole source is modeled at 
and made from two monopoles separated by a distance 
. Based on the regularized delta function 
chosen in (21), the separating distance is automatically set as 
.
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.
Note that two monopoles begin to overlap when 
. The corresponding superposition (blue line) will no longer satisfy the integration rule and cannot correctly model a dipole source.
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:The monopole source term 
in the wave equation (22) now becomes the summation of the two monopoles: 
.
The amplitude in the animation depends on the selection of the regularized delta function 
and should not be confused with the monopole source strength 
.
The analytical sound pressure amplitude of a dipole source is given by [23]:
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 was 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).
The amplitude in the animation depends on the selection of the regularized delta function 
and should not be confused with the dipole source strength 
.
The resulting sound field is not isotropic but is symmetric about the axis connecting the two monopoles.
Analytically, the dipole source 
is identical to the pair monopoles with a corresponding source strength. However, since the regularized delta function 
was applied 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 
.
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 equation is the combination of (28) and (29).
As an example, consider a 1D model of time-harmonic waves.
Note that when the model is built with 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 conditions, the specific acoustic impedance 
should be defined first. The specific acoustic 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 acoustic impedance can be shown to be equal to the product of the medium’s 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):
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.
A message is generated to indicate that the coupled first-order wave equation is convection dominated. This is as expected, since all the spatial derivative terms in (32) and (33) are first order only, representing the advection transport (or convective transport) [34] of the sound wave.
More information about this issue can be found in the "Finite Element Usage Tips" tutorial.
It is instructive to compare the result with the one computed with a second-order wave equation (35).
Note that the error, memory and time consumption are greatly reduced when the first-order system of PDEs is used. This is very helpful when solving large-scale acoustic models.
For a 2D version of the wave equation modeled as a set of coupled first-order equations, refer to here. Since the example presented there utilizes perfectly matched layers (PML), two extra factors 
and 
are introduced (known as PML absorption coefficients) in the governing equations (36). Once 
and 
are set to zero, the equation (37) will render into a coupled first-order wave equation (38) and (39) in the simulation domain, outside the PML region.
Sound Propagation in Lossy Media
To model sound propagation in lossy media such as porous materials, the wave equation (40) needs to be modified in order to accommodate for an attenuation of the sound waves [41]. 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 and 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 [42]. Note that equation (43) simplifies to the lossless wave equation (44) when 
, 
and 
.
Since the effective bulk modulus 
depends on a specific frequency 
, the lossy media model (45) is constrained to simulate single-frequency waves only, that is, a time-harmonic wave. To model a time-inharmonic wave that is composed of multiple frequencies, the wave should be decomposed by a Fourier transform before analysis.
As an example, 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.
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 (46) [47].
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 four types:
- Dirichlet-type boundary conditions. This type of boundary condition specifies the sound pressure
at the boundary and can be modeled with a DirichletCondition.
- Neumann-type boundary conditions. This type of boundary condition specifies the sound particle velocity
at the boundary and can be modeled with a NeumannValue.
- Robin-type boundary conditions. This type of boundary condition specifies the relation between time and normal derivatives of the sound pressure at the boundary and can modeled with a NeumannValue, since Robin-type boundary conditions are technically generalized Neumann boundary conditions.
- Periodic boundary conditions. This type of boundary condition maps the sound pressure
from one part of the boundary to another part and can be modeled with PeriodicBoundaryCondition.
Under these four 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, NeumannValue and PeriodicBoundaryCondition. 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 behavior of the boundary condition on the right-hand side.
Generally speaking, in solver algorithms, a NeumannValue[g-q p(t,X) ,X∈Γb] is used to specify the flux over the boundary 
such that 
holds.
However, in acoustic models, the dipole sources 
can only be specified within the domain and will always equal zero on any part of the boundary, and that leads to 
.
Impedance Boundary Conditions
Formulation
With a specified impedance 
on the boundary 
, the impedance boundary condition is given by:
Derivation
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 behavior, a tube filled with two sound media 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 pulse with standard deviation 
is used in the following example.
The region is split at 
into two media. The left medium 
is cold air 
and the right medium 
is hot air 
.
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.
The sound pressure field animation shows a sound wave moving from left to right, and two media 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 acoustic impedance. The specific acoustic 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 acoustic impedance can be shown to be equal to the product of medium’s density and the speed of sound [48]:
For a sound wave propagating from one medium to another, the amplitude relations can be expressed in specific acoustic impedance 
as:
where 
, 
and 
are the amplitude of the reflected, the incident and the transmitted wave, respectively. 
and 
are the specific acoustic impedance of the medium 1 and 2, respectively [49].
With these parameters, the solution can be verified.
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 (50) into the momentum conservation equation (51) 
, 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 (52). 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.
of medium 
, not including the part with medium 2 included in the 
used above:The animation shows the same result as the previous example, just with a smaller simulation domain.
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:
- For the case where both media are the same,
can be inserted in the equations and absorbing boundary conditions are obtained:
- For the case where the impedance at the boundary approaches infinity,
is inserted in the equations and sound hard boundary conditions are obtained:
- For the case where the impedance at the boundary approaches zero,
is multiplied on both sides of (53) and the limit is taken as 
approaches 0. Then the resulting equations are consistent with the definition of sound soft boundary conditions that 
. Note that sound soft boundary conditions should be directly modeled with DirichletCondition.
Absorbing Boundary Conditions
Formulation
With a specific type of incident wave and the distance 
between the wave origin to the boundary 
, the absorbing boundary condition is given by:
Derivation
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 makes 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 acoustic impedance 
of the domain. That is, an ABC is a special case of an impedance boundary condition that makes 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 acoustic impedance is 
. Substituting 
into the impedance boundary condition (54), then the absorbing boundary condition can be written as:
To accommodate for different kinds of waves, equation (55) can be generalized [56] with an additional term 
:
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, an infinitely long tube is shown, with the computational domain set from 
to 
. To model the continuation of the domain, an absorbing boundary condition is added at the right end.
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 how much an incoming wave is reflected back into the simulation domain. Ideally, nothing is reflected back.
The numerical error is due to the discretization of the finite element method. Note that the error did not 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.
is used:
See this note about improving the visual quality of the animation.
Note the reflection of the 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
Formulation
For a wall boundary 
, the sound hard boundary condition is given by:
Derivation
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 [57], the momentum conservation equation (58) relates the sound pressure 
with the sound particle velocity 
:
Combining (59) and (60), 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, a tube with one end closed is made use of.
at the closed end of the tube at 
: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 is the default boundary condition used if no boundary condition is specified at a given boundary.
This sound pressure field animation shows a sound wave moving from left to right toward 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
Formulation
With a specified sound particle velocity 
on the boundary 
, the normal velocity boundary condition is given by:
Derivation
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 
. This will ensure that the sound pressure will always be in phase with the particle velocity on a boundary.
Recall that the momentum conservation equation (61) relates the sound pressure 
with the sound particle velocity 
.
Inserting (62) into (63), 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.
The simulation begins with initial conditions that are at 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
.Note that this animation is quite similar to 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 acoustic impedance, so the perturbation on the sound pressure will also cause the perturbation on the sound particle velocity. When the specific acoustic impedance 
is known, a normal velocity boundary can be replaced by a pressure source boundary.
Sound Soft Boundary Conditions
Formulation
The sound soft boundary condition is given by:
Derivation
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 behavior, 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.
at the right open end of the tube: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
Formulation
With a specified sound pressure 
on the boundary 
, the pressure source boundary condition is given by:
Derivation
A pressure source boundary condition is spoken of 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.
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. An absorbing boundary condition is added at the left end to avoid wave reflection:
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 source on the boundary can be modeled by replacing the DirichletCondition with an equivalent monopole source on the right-hand side of the wave equation (64). Monopole sources are explained in the section entitled Monopole Sources.
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 acoustic impedance 
. Substituting 
into (65), the equation can be expressed in terms of the specified sound pressure 
as:
As an example, a time-dependent pressure source 
is modeled with the NeumannValue setting (66).
Note that the Dirichlet model is more accurate because Dirichlet conditions enforce 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 other Neumann value–based boundary conditions.
Radiation Boundary Conditions
Formulation
With a specified incident sound pressure 
and a wave direction vector 
on the boundary 
, the radiation boundary condition is given by:
The relation among the boundary normal vector 
, the wave direction vector 
and the wave incident angle 
is illustrated below:
Therefore, equation (67) can also be expressed with a wave incident angle 
as:
Note that when applying a radiation boundary in the 1D domain, the incident angle 
is always zero (i.e. normal incidence). Then equation (68) can be simplified to:
A 1D radiation boundary example is shown in the following section. A 2D example, which demonstrated a radiated sound wave with oblique incidence, is presented in the appendix section: Radiated Sound Waves with Oblique Incidence.
Derivation
In this section, the radiation boundary condition is derived for the 1D domain. For a complete derivation in the 2D and 3D domain, refer to the study from Givoli and Neta [69].
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 (70) and (71) and inserting (72) gives:
Recall that on the boundary, the amplitude of the incoming wave is specified by 
. Substituting this relation into the right-hand side of (73) yields the radiation boundary condition:
To illustrate this behavior, 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.
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 moves neither right nor left but simply oscillates in time.
This type of wave is called a standing wave and will be explained in more detail in the separate tutorial entitled "Acoustics in the Frequency Domain".
Periodic Boundary Conditions
Formulation
Given a function 
that maps the sound pressure 
from the source boundary 
to the periodic boundary 
, the periodic boundary condition can be written as:
Derivation
A periodic boundary condition is used to model an acoustic wave in a spatially periodic domain. That is, the information obtained from one boundary (referred to as the "source boundary 
") can be mapped to another boundary (referred to as the "periodic boundary 
") with a relating function 
. The periodic boundary condition is set by the PeriodicBoundaryCondition in the acoustic PDE model.
As an example, an acoustic wave propagating in a circular tube is shown. With the usage of the periodic boundary condition, it is possible to perform the simulation in a 1D domain.
The circular tube is converted into a 1D model of the length 
equal to the perimeter of the tube. The periodic boundary 
is set on the left end so that when the sound wave passes through the right side of the domain 
, it reappears on the left side with the same magnitude.
In this example, an initial acoustic wave is set with the error function Erf. Note that the derivative of the initial condition is not zero but chosen such that the initial wave travels to the right.
On the right end an absorbing boundary condition is added to avoid the wave reflection.
The animation shows a right-traveling sound wave within a spatially periodic domain. As the wave passes through the right boundary, it reappears on the left side because of the periodicity of the domain set up with the periodic boundary condition. Since the sound medium is assumed to be acoustic lossless, the magnitude of the wave has been kept at a constant level at all times.
Perfectly Matched Layer
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 [74]. 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 PMLs can be found in the appendix section entitled Perfectly 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 equation with first-order-only spatial derivatives. Then the entire process of generating a PML can be performed by a single transformation (75) 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 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, a simulation of a 1D domain is considered from 
to 
using a computational domain that reaches from 
to 
and a PML at the right from 
to 
.
Take air as the sound medium, and for an one-dimensional plane wave the specific acoustic impedance 
is 
[76].
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.
As seen, 
in the original non-PML region.
A message is generated to indicate that the coupled first-order wave equation is convection dominated. This is as expected, since all the spatial derivative terms in (77) are first order only, representing the advection transport (or convective transport) [78] of the sound wave. This warning message can be switched off by evaluating Off[NDSolveValue::femcscd].
More information about this issue can be found in the "Finite Element Usage Tips" tutorial.
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.
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 
.
To implement the PML on a 2D wave equation, two auxiliary field variables 
and 
are introduced [79]. The equation becomes:
Here 
and 
are sound particle velocity components in space.
is used:
See this note about improving the visual quality of the animation.
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 [80]. The key property of a PML that distinguishes it from ordinary absorbing material is that it is designed so that waves propagated 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 
-direction is expressed as:
The concept of a PML is to perform a coordinate transformation so that all the 
terms in wave expression (81) 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, a simulation of a 1D domain from 
to 
is considered using a computational domain that stretches from 
to 
with a PML at the right from 
to 
. To achieve this, a linearly growing imaginary part 
with the slope 
is added in the PML.
With the PML coordinate transformation (82), the wave expression (83) transforms to:
Since there is an extra term 
in the new wave expression (84), the wave with positive wave number 
will exponentially decay in the absorbing layer where 
.
Note that the wave is unchanged outside the PML (
) where 
.
for the artificial imaginary part 
:Now the wave is attenuated rapidly and decays to nearly zero before 
. In other words, the greater attenuation rate (larger 
value) allows a thinner PML to be used.
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 First-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 (85) are replaced by the complex-valued variable 
defined as:
In order to make the PML attenuation rate independent of the wave frequency 
[86], 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:
However, the original wave equation (87) contains the second-order spatial derivative. So it is more convenient to split the wave equation into coupled differential equations with only the first-order spatial derivative, then the entire process of PML can be performed by a single transformation (88) 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 perform the PML coordinate transformation (89) on (90) 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:
(91) can be written as a system of equations, and in fact, this transformation is done internally in NDSolve. Introduce the following auxiliary equation for the transformation:
Substitute (92) into (93) and use (94) to obtain an equivalent system of first-order equations for the original wave equation (95):
Next, the substitution from (96) is done and the initial conditions are modified.
The important point is to note how the time derivative in the Neumann value changed by applying (97) to the auxiliary dependent variable. In this case, the time derivative 
in the Neumann value is replaced by a Neumann value 
on the auxiliary variable.
Radiated Sound Waves with Oblique Incidence
To model radiating sound waves with an arbitrary incident angle 
, the general equation for the radiation boundary condition can be used:
Here 
denotes the incident sound pressure on the boundary 
.
As a demonstration, consider the following 2D example. A time-harmonic sound wave is radiating from the bottom boundary 
with an incident angle of 
. The wave then travels through the domain 
and leaves from the top boundary 
.
For an irregular geometry, it is convenient to make use of element markers for boundary condition specification. More information on markers and their generation in meshes can be found in a separate tutorial: "Element Mesh Generation: Markers".
See this note about improving the visual quality of the animation.
The animation shows a non-normal sound wave radiated from the bottom boundary 
at 
. The wave then travels through the domain 
and leaves from the top boundary 
.
Nomenclature
The following table summarizes symbol names used and their meaning in this tutorial.
| Symbol | Description | Unit |
| ρ | density of a medium | [kg/m3] |
| c | speed of sound in a medium | [m/s] |
| p | sound pressure | [Pa] |
| pmax | maximum value of sound pressure | [Pa] |
 | sound pressure amplitude function | [Pa] |
 | specified boundary pressure | [Pa] |
| t | time | [s] |
| tend | simulation end time | [s] |
| X | position vector | [m] |
| s | direction switch | N/A |
| F | optional dipole source | [N/m3] |
 | dipole source strength | [N/m3] |
| Q | optional monopole source | [1/s2] |
 | monopole source strength | [1/s2] |
| d | separation distance of dipole source | [m] |
| λ | wavelength of sound | [m] |
| Ω | simulation domain | [m] |
| k | wave number | [rad/m] |
| f | sound wave frequency | [Hz] |
| ω | sound wave angular frequency | [rad/s] |
| δ | Dirac delta function | N/A |
 | regularized delta function | N/A |
| γ | regularization parameter | [m] |
| h | mesh spacing | [m] |
| Xs | sound source location | [m] |
| α | attenuation factor | [m2/(s·N)] |
| ϕ | porosity | N/A |
| VV | void volume | [m3] |
| VT | total volume | [m3] |
| Rf | flow resistivity | [kg/(m3·s)] |
| β | standard deviation of a Gaussian pulse | [m] |
| ζ | sound particle displacement | [m] |
| v | sound particle velocity | [m/s] |
 | specified boundary velocity | [m/s] |
| T | temperature | [K] |
| Z | specific impedance | [Pa·s/m] |
| Zb | boundary impedance | [Pa·s/m] |
| Ar | amplitude of reflected wave | [Pa] |
| Ai | amplitude of incident wave | [Pa] |
| At | amplitude of transmitted wave | [Pa] |
 | complex spatial variable of PML | [m] |
| fp | PML model parameter | [m] |
| mf | slope of PML artificial imaginary part | N/A |
| σ | absorbtion coefficient of PML | [rad/(s·m)] |
| σmax | maximum value of absorbtion coefficient | [rad/(s·m)] |
| ψ | arbitrary function | N/A |
References
1. Nolte, B. and Marburg, S. Computational Acoustics of Noise Propagation in Fluids: Finite and Boundary Element Methods. Springer-Verlag, 2008.
2. Heutschi, K. Lecture Notes on Acoustics I. Swiss Federal Institute of Technology Zurich, 2016.
3. Russell, D. "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, K. "Finite Element Analysis of 2D Acoustics." Kai’s Thoughts, 2013. https://blog.kaistale.com/?p=574
5. Fahy, F. Foundations of Engineering Acoustics. Academic Press, 2001.
6. Johnson, S. "Notes on Perfectly Matched Layers (PMLs)." MIT, 2010.
7. Bilbao, S. and Hamilton, B. "Directional Source Modeling In Wave-Based Room Acoustics Simulation." IEEE, 2017.
8. Peskin, C. "The Immersed Boundary Method." Cambridge University, 2002.
9. Russell, D., Titlow, J. and Bemmen, Y.-J. "Acoustic Monopoles, Dipoles and Quadropoles: An Experiment Revisited." American Journal of Physics, 67, 660, 1999.
10. Zeleny, E. "Acoustic Multipoles," Wolfram Demonstrations Project, 2013. http://demonstrations.wolfram.com/AcousticMultipoles/
11. Zeleny, E. "Longitudinal and Transverse Waves," Wolfram Demonstrations Project. http://demonstrations.wolfram.com/LongitudinalAndTransverseWaves/
12. Moore, G. "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, F. "The Medium-Frequency Range in Computational Acoustics: Practical and Numerical Aspects." Journal of Computational Acoustics, Vol.11, No. 2 175-193, 2003.
15. Vita, M. "The Wave Equation with a Source." Oklahoma State University.
16. Maderuelo-Sanz, R. 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. Cox, T. and D'Antonio, P. Acoustic Absorbers and Diffusers: Theory, Design and Application. Spon Press, 2004.
19. Engquist, B. and Majda, A. "Absorbing Boundary Conditions for the Numerical Simulation of Waves." Mathematics of Computation, Vol. 31, No. 139 629-651, 1977.
20. Duran, D. Numerical Methods for Fluid Dynamics - With Applications to Geophysics, 2nd ed., Springer, 2010.
21. Givoli, D. and Neta, B. "High-Order Non-reflecting Boundary Scheme for Time-Dependent Waves." Journal of Computational Physics, Vol. 186, 24–46, 2004.