Acoustics in the Time Domain
|Acoustic Boundary Conditions||References|
|Perfectly Matched Layer (PML)|
Acoustics is the field of physics that models sound by changes in pressure. The changes in pressure are described by a wave equation. This tutorial gives an introduction to modeling sound with the wave equation in the time domain and presents various aspects of the modeling process. The modeling process results in partial differential equation (PDE) models that are solved with NDSolve. Furthermore, in this tutorial different types of sound sources are introduced as well as an overview of how various real world sound barriers can be modeled with the available PDE boundary conditions.
In this tutorial time domain modeling is used to investigate how sound waves behave in time. It should be noted that sound can also be modeled in the frequency domain. Employing the frequency domain for modeling illustrates how sound waves are distributed over a range of frequencies. Frequency domain modeling is widely used when analyzing acoustic systems with a frequency dependency, such as noise filters.
Details about the modeling of sound in the frequency domain can be found in a separate tutorial entitled Acoustics in the Frequency Domain.
The symbols and corresponding units used throughout this tutorial are summarized in the Nomenclature section.
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 . The following animation illustrates particle motions of longitudinal waves and transverse waves .
The animation illustrates how the highlighted particle oscillates horizontally on longitudinal waves, but fluctuates vertically on transverse waves.
Defining a function that represents the spatial terms of an acoustic model will make the setup for the wave equation more convenient.
Various sections in the documentation explain the use of inactive PDE operators. Please refer to "Numerical Solution of Partial Differential Equations".
To make use of specific material parameters in the equation we extract relevant data from the ThermodynamicData.
Sound waves can be classified into two categories based on their form of time dependence. Time harmonic waves have sinusoidal time variation with a single frequency. Time inharmonic waves have arbitrary time variation with multiple frequencies.
In acoustics simulations the wavelength of a sound wave needs to be resolved by a sufficiently fine mesh in order to get an accurate numerical solution of the governing partial differential equation. Here we create a function to set the max cell measure for both time harmonic and time inharmonic waves.
Throughout this tutorial we will use a particular example frequency for time harmonic waves and a Gaussian pulse as an example of time inharmonic waves.
The default resolution of the mesh is set to 20 elements per wavelength , which means the max cell measure would be twenty times smaller than the sound wavelength . Typically this is sufficient to resolve the waves accurately . 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 .
Define a function named AcousticOptions that gives NDSolve options appropriate for solving acoustic problems.
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.
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).
An arbitrary wave can be decomposed into an infinite sum of harmonic waves with discrete frequencies . 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.
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.
Use Manipulate to investigate a wave with adjustable frequency 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 increase when selecting a higher frequency . The number of phase changes per unit distance increase as the wave number increase. The wave propagation can be observed by playing the time slider.
When is substituted in the source free wave equation (10), the equation simplifies to an inhomogeneous Helmholtz equation that is a second order time independent differential equation.
In the case where there are no sources this PDE system becomes an eigenvalue problem that can be solved with NDEigensystem. The set of eigenvalues then determines the natural frequencies of the corresponding acoustic system.
The details of the sound modeling with the Helmholtz equation can be found in a separate tutorial entitled Acoustics in the Frequency Domain.
Next, we build a 1D time harmonic wave model using the wave equation (11) and solve it numerically using NDSolveValue.
A sinusoidal source on the left end of will produce a time harmonic wave. The source can be added using DirichletCondition.
A point is shown in red to demonstrate that the amplitude (pressure) of the time harmonic wave at a particular value of oscillates sinusoidally with the same frequency as the incoming wave.
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 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 . 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
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.
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.
Determine appropriate mesh spacing using the function introduced in Model Parameter Setup.
Define the monopole source as a function of , , and the location . Here we use the mesh spacing as the regularization parameter .
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 .
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.
Use as the regularization parameter to ensure that the two monopoles do not overlap. The negative sign of indicates that the two monopoles are in the opposite phase.
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 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 .
Show the difference between the dipole model and pair monopoles model (blue) superimposed on the wave scaled to fit on the same vertical scale (light gray).
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 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 :
Note that when we build the model with the coupled first order wave equations, extra initial and boundary conditions are required for the sound particle velocity
To set up the extra initial and boundary condition, the specific impedance should be defined first. The specific impedance denotes the ratio of the sound pressure to the sound particle velocity as the wave moves through the medium, which is written as:
For a plane wave the specific impedance can be shown to be equal to the product of the medium’s density and the speed of sound . A detailed explanation of impedance can be found in the section entitled Impedance Boundary Conditions.
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.
It is instructive to compare the result with the one computed with a second order wave equation (32).
Solve the second order PDE and monitor the time/memory usage.
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.
To model sound propagation in lossy media such as porous materials, the wave equation (33) needs to be modified in order to accommodate for an attenuation of the sound waves . 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 . Note that equation (36) simplifies to the lossless wave equation (37) when , , and .
Since the effective bulk modulus depends on a specific frequency , the lossy media model (38) is constrained to simulate single frequency waves only, that is, a time-harmonic wave. To model a time-inharmonic wave that composes of multiple frequencies, the wave should be decomposed by a Fourier transform before analysis.
As an example we consider a harmonic sound wave at frequency , which propagates through a tube filled with a porous absorber. The porosity and flow resistivity of this material are given by and , respectively.
An absorbing boundary condition is added on the right to avoid a reflection of the wave. Note that the NeumannValue is set to to accommodate for the modified wave equation (39) .
A harmonic sound source at frequency is added on the left end using DirichletCondition.
For comparison set up an analytical time harmonic wave propagating in a lossless media.
In the lossy media the sound wave propagates in a lower speed while being attenuated within the domain.
The most common boundary conditions in acoustics can be modeled with DirichletCondition and NeumannValue, and can be categorized in the following three 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.
The following section describes several physical boundaries common in acoustics and how they can be modeled with the use of DirichletCondition and NeumannValue. For this purpose the boundary condition currently discussed is always on the right hand side of the simulation domain. In some examples additional boundary conditions are on the left hand side to better demonstrate the behaviour of the boundary condition on the right hand side.
When a sound wave propagates from one medium into another, part of the sound wave will be reflected at the interface and part will be transmitted across the interface. To illustrate this behaviour a tube filled with two sound mediums is considered in the next example.
For a tube aligned with the axis and a diameter much smaller than the sound wave length , it is then reasonable to assume that all points in the , plane have the same sound pressure value and the same sound particle velocity . This wave type is called plane wave and can be treated as a one-dimensional wave.
To better illustrate an impedance boundary condition a Gaussian pulses with standard deviation is used in this example.
The region is split at into two media. The left medium is cold air and the right medium is hot air .
Specify speeds of sound and densities with ThermodynamicData.
In 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 impedance. The specific impedance is the ratio of the sound pressure to the sound particle velocity as the wave moves through the medium, which is defined as:
For a plane wave the specific impedance can be shown to be equal to the product of medium’s density and the speed of sound .
For a sound wave propagating from one medium to another, the amplitude relations can be expressed in specific impedance as
where , and are the amplitude of the reflected, the incident and the transmitted wave, respectively. and are the specific impedance of the medium 1 and 2, respectively .
If one is not interested in the solution in, say, medium then there is no need to include that region in the solution. The medium interface can be modeled by specifying the impedance of medium at the boundary at . This is called an impedance boundary condition. An impedance boundary condition can be realized with a generalized Neumann boundary condition which is also known as a Robin boundary condition.
Substituting the definition of impedance (43) into the momentum conservation equation (44) , then the impedance boundary condition can be expressed in the boundary impedance as
Rebuild the example with known boundary impedance , so the NeumannValue is set to . Note that the coefficient in needs to match the NeumannValue specification . Here and are 0 and matches and this is then consistent with the wave equation (45) we use. Also, the medium boundary is at the right end where . The meaning of time derivatives in NeumannValue is explained in more detail in the section entitled Neumann Value on Time Derivative.
Set up the impedance boundary condition with NeumannValue.
Solve the sound pressure field only in the region of medium , not including the part with medium 2 included in the used above.
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.
- 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 (46) 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.
Typically a simulation domain that extends to infinity is not a feasible option in a simulation. Absorbing boundary conditions are a methodology used to model boundaries at infinity. An absorbing boundary condition (ABC) works by absorbing an incoming wave and thus make the model behave as if it had infinite extent. ABCs are not the only way to model simulation domains with infinite extent. As an alternative approach to ABCs perfectly matched layers (PML) may be used.
Since an absorbing boundary is a continuation of the computational domain, the boundary impedance is equal to the specific impedance of the domain. That is, an ABC is a special case of an impedance boundary condition that make use of the time derivative of the dependent variable within the NeumannValue. The meaning of time derivatives in NeumannValue is explained in more detail in the section entitled Neumann Value on Time Derivative.
Recall that for a plane wave the specific impedance is . Substituting into the impedance boundary condition (47), then the absorbing boundary condition can be written as
To accommodate for different kinds of waves, equation (48) can be generalized with an additional term : 
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.
Set the initial conditions of the sound pressure field. A 2D Gaussian pulse with the standard deviation is used.
Solve the PDE with NDSolveValue and monitor time/memory usage.
To visualize the wave reflection, we set up a legend bar and ContourPlot options scaling from to .
Visualize the truncated solution with the vertical PlotRange set sufficiently small to see the reflections.
Note that the reflected wave at the four corners where the incident angles are not normal to boundaries. With a different geometry the reflections may be somewhat smaller.
To reduce the wave reflection in 2D or 3D, a perfectly matched layer (PML) can be applied as an alternative method to model an absorbing boundary. A PML is designed such that waves that impact boundaries in a non-normal way are absorbed within it. Further details on PML can be found in the section entitled Perfectly Matched Layer (PML).
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 , the momentum conservation equation (51) relates the sound pressure with the sound particle velocity .
Combining (52) and (53), it is found that a sound hard boundary condition also implies that the partial derivative of the sound pressure in direction normal to the boundary is zero.
To model a sound hard boundary condition the NeumannValue is set to 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 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.
A boundary where the sound particle velocity normal to the boundary is specified and not equal to zero is called a normal velocity boundary. A familiar example of a normal velocity boundary is a vibrating mechanical structure that produces sound waves into the surrounding medium.
The negative sign is added in front of to indicate that the sound particle velocity is specified opposite to the outward normal . By doing so we ensure that the sound pressure will always be in phase with the particle velocity on a boundary.
Recall that the momentum conservation equation (54) relates the sound pressure with the sound particle velocity .
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.
Set up the normal velocity condition at the right end with NeumannValue.
The simulation begins with initial conditions are rest and the normal velocity boundary condition generates a sinusoidal wave on the right hand boundary that propagates to the left. The small wiggles ahead of the wave are artifacts due to the spatial discretization. This effect can be reduced by setting a higher mesh resolution in AcousticOptions.
Note that this animation is quite similar with the example shown in Pressure Source Boundary Conditions. That is because the sound pressure and particle velocity have linear dependence with each other given by , where is the specific impedance, so the perturbation on the sound pressure will also cause the perturbation on the sound particle velocity. When the specific impedance is known a normal velocity boundary can be replaced by a pressure source boundary.
On a sound soft boundary the medium pressure at the boundary is set equal to an ambient reference pressure such that . This type of boundary is also called sound free boundary since the sound particles are free to move at the boundary, which means that the sound particle velocity is not required to be zero. The momentum conservation equation relates the sound particle displacement and sound pressure . While , is not necessarily equal to zero.
At a sound soft boundary the impedance since the sound pressure . That is, the sound soft boundary condition is a special case of impedance boundary conditions when the boundary impedance is set to zero:
To illustrate this behaviour a tube with one end open is considered in the next example. The open-ended side is treated as a sound soft boundary since there is no constraint to limit the sound particles movement.
To model the sound soft boundary condition a DirichletCondition is set to at the right open end of the tube.
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 .
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.
Set up the pressure source boundary with DirichletCondition.
The simulation begins with an undisturbed domain and the transient pressure variation generates a sinusoidal wave on the right hand boundary that propagates to the left.
Note that a pressure source generates waves normal to the boundary while a monopole source radiates waves isotropically. So for 1D cases, a pressure sources on the boundary can be modeled by replacing the DirichletCondition with a equivalent monopole source on the right hand side of the wave equation (57). Monopole sources are explained in the section entitled Monopole Sources.
A pressure source can also be modeled with a NeumannValue. Doing so has the advantage that on the same portion of the boundary the pressure source can be combined with other Neumann values based boundary condition models. In other words hybrid boundary conditions can be created by combining several Neumann values on the same part of the boundary. In this way, for example, a pressure source can be combined with an absorbing boundary to form a radiation boundary condition.
The NeumannValue setting for a pressure source can be derived from the normal velocity boundary conditions:
Recall that for a plane wave the specific impedance . Substituting into (58) the equation can be expressed in terms of the specified sound pressure as
As an example we model a pressure source with the NeumannValue setting (59).
Set up the pressure source boundary with NeumannValue.
Note that the Dirichlet mode is more accurate because Dirichlet conditions enforce a specific values at the boundary nodes. This means that modeling a pressure source with a Neumann value should only be done if there is a need to combine that pressure source with a other Neumann value based boundary conditions.
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.
Recall that on the boundary the amplitude of incoming wave is specified by . Substituting this relation into the right hand side of (63) yields the radiation boundary condition
To illustrate this behaviour a semi-infinite tube with one end closed is considered in the next example. A harmonic sound wave enters the computational domain from the right end. The wave then travels through the tube. The left end is then treated as a sound hard boundary. When the wave comes back to the right end it is absorbed. This process of being a pressure source and absorbing the reflected wave is the radiation boundary condition.
Set up the radiation boundary with NeumannValue at the right end. By default a sound hard boundary condition is implicitly used at the left end.
Solve the PDE with NDSolveValue. To show the effect of the radiation boundary the simulation time is extended a little.
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.
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 . The key property of a PML that distinguishes it from an ordinary absorbing boundary condition is that it is designed so that waves propagate from a non-PML region to a PML and do not reflect at the interface in any incident angle. That is, a PML acts like a non-reflecting absorbing material. Further details on PML can be found in the appendix section entitled Perfect Matched Layer Derivation.
To implement a PML two things need to happen. First the simulation domain needs to be enlarged. This extension is the region in which the PML is active. Second, a coordinate transformation of the PDE is done.
It is better to split the wave equation into a coupled differential equations with first-order only spatial derivatives. Then the entire process of generating a PML can be performed by a single transformation (65) of the original wave equation.
The coupled first-order wave equation after the PML coordinate transformation is
The absorption coefficient is a parameter to be chosen, and is set to increase linearly within the PML layer from to . A thinner PML requires a greater value of but large tends to increase numerical reflections.
A PML interface is only reflectionless when we are solving the exact, continuous wave equation. NDSolve, however, is a numerical method that discretizes the simulation domain and solves an approximation to the wave equation instead. In this numerically approximate case the wave will show reflections at the PML interface. To reduce the amount of reflection a finer grid or a smaller attenuation rate (smaller ) should be used.
For an example we consider a simulation of a 1D domain from to using a computational domain that from to with a PML at the right from to .
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.
Solve the PDE with NDSolveValue. To show the effect of PML the simulation time is extended a little.
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 . The equation becomes:
Set the initial conditions of the sound pressure field. A 2D Gaussian pulse with the standard deviation is used.
Solve the PDE with NDSolveValue and monitor time/memory usage.
To visualize the wave reflection, we set up the legend bar and the ContourPlot options scaling from to .
Visualize the truncated solution in the non-PML region with the PlotRange set to .
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.
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 . 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.
The concept of a PML is to perform a coordinate transformation that all the terms in wave expression (69) are replaced by a complex variable .
where the original serves as the real part of , and is an artificial imaginary part in an absorbing layer.
This way a wave will decay exponentially in PML but remain unchanged otherwise. The PML attenuation rate is controlled by the slope of the artificial imaginary part . For a thinner PML a greater attenuation rate is required, that is, a larger value has to be assigned. However, using a very large value can increase the numerical reflection of the PML, so there is a tradeoff between the PML width and 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.
Since there is an extra term in the new wave expression (72), the wave with positive wave number will exponentially decay in the absorbing layer where .
To increase the attenuation rate in PML, we can pick a greater slope value 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 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 (73) are replaced by the complex valued variable defined as:
In order to make the PML attenuation rate independent of the wave frequency , 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
However, the original wave equation (75) contains the second order spatial derivative. So it is more convenient to split the wave equation into coupled differential equations with only first-order spatial derivative, then the entire process of PML can be performed by a single transformation (76) of the original wave equations.
Recall that the coupled first-order wave equation is shown as
Now we perform the PML coordinate transformation (77) on (78), and apply the harmonic wave relation and . The equations become
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
An absorbing boundary condition can be modeled with a NeumannValue on the time derivative. This is best understood when looking at a simplified wave equation.
(79) can be written as a system of equations; and in fact this transformation is done internally in NDSolve. For the transformation we introduce the following auxiliary equation:
If we now substitute (80) into (81) and use (82) we obtain an equivalent system of first order equations for the original wave equation (83).
The important point is to note how the time derivative in the Neumann value changed by applying (85) to the auxiliary dependent variable. In this case the time derivative in the Neumann value is replaced by to become a Neumann value on the auxiliary variable.
|ρ||density of a medium||[kg/m3]|
|c||speed of sound in a medium||[m/s]|
|pmax||maximum value of sound pressure||[Pa]|
|sound pressure amplitude function||[Pa]|
|specified boundary pressure||[Pa]|
|tend||simulation end time||[s]|
|F||optional dipole source||[N/m3]|
|dipole source strength||[N/m3]|
|Q||optional monopole source||[1/s2]|
|monopole source strength||[1/s2]|
|Null||separation distance of dipole source||[m]|
|λ||wavelength of sound||[m]|
|Null||sound wave frequency||[Hz]|
|ω||sound wave angular frequency||[rad/s]|
|δ||Dirac delta function||N/A|
|regularized delta function||N/A|
|Null||sound source location||[m]|
|β||standard deviation of a Gaussian pulse||[m]|
|ζ||sound particle displacement||[m]|
|Null||sound particle velocity||[m/s]|
|Null||specified boundary velocity||[m/s]|
|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)]|