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.

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 the 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. 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 max edge length is set to 12 nodes per wavelength , which means that there will be at least twelve 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 we will typically be using Gaussian pulses with standard deviation .

In many examples we will be using an absorbing boundary condition (ABC) 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.

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:

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.

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 the time slider.

Next, we build a 1D time inharmonic wave model using the wave equation (13) and solve it 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 .

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 .

Modifying will change the width of the regularized delta function, however, for the spatial integral is always 1.

Next, we build a 1D time inharmonic wave model using the wave equation (19) and solve it numerically using NDSolveValue.

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.

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 .

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 loner satisfy the integration rule and can not correctly model a dipole source.

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 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).

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 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 .

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 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 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 we use the first order system of PDEs. This is very helpful when we are solving a large scale acoustic models.

For a 2D version of the wave equation modeled as a set of coupled first order equations please 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, 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 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.

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:

• Robin Type

• Impedance Boundary Conditions

• Absorbing Boundary Conditions

• Radiation Boundary Conditions

• Neumann Type

• Sound Hard Boundary Conditions - Walls

• Normal Velocity Boundary Conditions

• Pressure Source Boundary Conditions (Neumann model)

• Dirichlet Type

• Sound Soft Boundary Conditions

• Pressure Source Boundary Conditions (Dirichlet model)

• Periodic Type

• Periodic Boundary Conditions

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 behaviour 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 to zero on any part of the boundary and that leads to .

Impedance Boundary Conditions

With a specified impedance on the boundary , the impedance boundary condition is given by:

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 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 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 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) 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.

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 we insert in the equations and we obtain absorbing boundary conditions:

• For the case where the impedance at the boundary approaches infinity, we insert in the equations and we obtain sound hard boundary conditions:

• For the case where the impedance at the boundary approaches zero, we multiply on both sides of (53) and take the limit as approaches 0. Then the equations we obtain 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

With a specific type of incident wave and the distance between the wave origin to the boundary , the absorbing boundary condition is given by:

• Plane wave:

• Cylindrical wave:

• Spherical wave:

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 acoustic 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 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 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.

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.

The numerical error is due to the discretization of the finite element method. Note that the error didn’t 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.

See this note about improving the visual quality of the animation.

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 layers (PML).

Sound Hard Boundary Conditions - Walls

For a wall boundary , the sound hard boundary condition is given by:

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 we look at a tube with one end closed.

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.

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

With a specified sound particle velocity on the boundary , the normal velocity boundary condition is given by:

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 (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 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 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

The sound soft boundary condition is given by:

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.

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

With a specified sound pressure on the boundary , the pressure source boundary condition is given by:

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.

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.

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 (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 we model a time dependent pressure source with the NeumannValue setting (66).

Note that the Dirichlet mode 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

With a specified incident sound pressure and a wave direction vector on the boundary , the radiation boundary condition is given by:

The relation between 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 the equation (68) can be simplified to:

A 1D radiation boundary example will be 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.

In this section we will derive the radiation boundary condition for the 1D domain. For a complete derivation in 2D and 3D domain please 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 incoming wave is specified by . Substituting this relation into the right hand side of (73) 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.

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.

Periodic Boundary Conditions

Given a function that maps the sound pressure from the source boundary to the periodic boundary , the periodic boundary condition can be written as:

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 as the "source boundary ") can be mapped to another boundary (referred 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 we look at an acoustic wave propagating in a circular tube. 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 equals to the perimeter of the tube. The periodic boundary is set on the left end so that when the sound wave pass through the right side of the domain , it re-appears 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 pass through the right boundary it re-appears 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 (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 [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 PML 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 equations 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 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 .

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.

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 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.

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 (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 it’s 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.

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 .

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.

Note that there is no reflection at the PML interface (), so it acts like a non-reflecting absorbing material.

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], 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 (87) 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 (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 we 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. For the transformation we introduce the following auxiliary equation:

If we now substitute (92) into (93) and use (94) we obtain an equivalent system of first order equations for the original wave equation (95):

Next, we perform the substitution from (96) and modify the initial conditions.

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 symbols 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, 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. Kai’s 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.

20.  D. Duran. Numerical Methods for Fluid Dynamics - With Applications to Geophysics, second edition, Springer, 2010.

21.  D. Givoli and B. Neta. High-order Non-reflecting Boundary Scheme for Time-dependent Waves. Journal of Computational Physics, Vol. 186, 24–46, 2004.

Related Tech Notes

▪
• Acoustics in the Time Domain
▪
• Acoustics in the Frequency Domain
▪
• Acoustics Model Collection
▪
• PDEModels Overview