Introduces technical terms related to sound and audio technology, essential physics in sound and audio research and development, terms used in mathematics, and terms on technology that are attracting attention as a trend of sound and audio technology.
Keywords: Acoustic analysis, finite difference time area law, FDTD
In recent years, there are three main types of methods used for acoustic numerical analysis. The Finite Element Method (FEM), the Boundary Element Method and the Finite Difference Time Domain Method (FDTD) are the Finite Difference Time Domain Method. The previous two are often used in the prediction of steady state (Note 1), while the FDTD method is a method of calculating all times from the initial state to the specified time or steady state. There is.
First, I will briefly introduce how the FEM method and the BEM method are applied to speaker development. Since the FEM method can efficiently analyze solids and fluids, it can verify how the sound generated by the combination of the shape of the diaphragm and the material changes (see Fig. 1 (A). ) Since the BEM method can calculate the acoustic characteristics on the boundary, the orientation frequency characteristics of the speaker can be verified by calculating the acoustic characteristics of 1m or 2m, mainly speakers (Fig. 1 (B)). )reference).
Figure 1 FEM method and BEM method used for speaker development
The FDTD method is still rarely used for speaker development, and is widely used in the architectural sound field that emphasizes sound wave propagation [1, 2].
The FDTD method is a method that is originally used to analyze electromagnetic waves. Based on the Maxwell equation, the time change of electromagnetic distribution in the observation space can be calculated using a relevant formula between magnetic fields and electric fields. The same approach allows the sound analysis (Note 2) by the FDTD method for the same wave phenomenon.
Note 1: Standing state does not change depending on time.
Note 2: Acouter analysis calculates the sound wave shape at arbitrary points in the space, and also requires each sound indicator (frequency characteristics, reverberation time, sound insulation performance, etc.).
The interesting thing about the sound field analysis by the FDTD method is that it can seek phenomena such as reflection, diffraction, and interference when the sound waves arrive in an obstacle. By visualizing this time change as a video, you can clearly indicate the state of sound wave propagation.
Here is an example of sound wave propagation analysis. Figure 2 and 3 assume a square meter spatial area, and set the area around it, such as a free sound field, is set. In Figure 2, when an impulse sound source is given to the center of the space, it is possible to confirm that the sound wave is propagated as a circle centered on the sound source. Further, in Fig. 3, a square obstacle is placed in the lower right to observe the phenomenon of reflection and diffraction by obstacles.
Figure 2 Acoustic analysis video by FDTD method
Figure 3 音 Acoustic analysis video by the FDTD method: If there is an obstacle in the lower right
Next, I will introduce the basic elements required for acoustic analysis using the FDTD method. It is the setting of the formula, stable conditions, boundary conditions and sound sources. There are more complicated methods for these elements, but all of them are applied to calculus and calculator theory, so I will omit it in this article.
(1) Formatization of the mathematical model of sound wave propagation and the FDTD method First, I will explain from the mechanism of sound wave propagation. When vocal cords and vibration plates of speakers begin to vibrate as the source of sound, the surrounding air particles begin to move due to this vibration. With this movement, the air density changes and pressure (sound pressure) changes (state 1). And due to pressure changes, the nearby air particles also begin to move (state 2). The sound waves are formed by the repetitive periodic movement of this air particles, and the sound waves are propagated from the sound source.
Express each of the above conditions with a simple one -dimensional mathematical model.
Condition 1: As shown in Fig. 4. A, the air particles are densely $ $ $ and speed $ U $ on the air grid to be observed, and the speed of $ U+du $ leaks at this instant lattice. The mass change rate can be expressed as follows. \ begin {gather} \ Frac {d \ rho^\ prime} {dt} dx = \ rho_0 u -\ rho_0 (u + du) = - \ rho_0 du \ nonumber \\ [10pt] \ Frac {d \ rho^\ prime} {dt} = - \ rho_0 {du} {dx} \ end {gather} Here $ ρ '$ is the air density in the lattice. In addition, since there is a relationship between $ p = ρ'c^2 $ from the air condition equation, the equation (1) is (2). \ begin {gather} \ Frac {dp} {dt} = - \ rho_0 c^2 {du} {dx} \ end {gather} Here $ C $ is the speed of sound in the air.
Condition 2: If a pressure difference $ DP $ occurs on both sides of the air grid as shown in Fig. 4.B, the air particles in the lattice move at $ U $ in the second law of the exercise. The relationship between the pressure and speed is as follows. \ begin {gather} \ Frac {dp} {dx} = -\ rho_0 {du} {dt} {dt} \ end {gather}
Figure 4 Equal equations and exercise equations, one -dimensional models
Equation (2) and (3) are called continuous equations of sound wave propagation and exercise equation. The FDTD method can convert these two differential equations to differential equations and calculate sound wave propagation only in four arithmetic operations. The difference formula (2) is as follows. \ begin {gather} \ Frac {p_i^{n+1} -p_i^n} {\ delta t} = - \ rho_0 c^2 {u_ {i+1/2}^{n+1/2} -U_ {i {i {i_ {I_U_ {{{{{{{{{{{{{} -1/2}^{n+1/2}} {\ delta x} \ nonumber \\ [12pt] P_i^{n+1} = p_i^n - \ delta t \ Left ( \ rho_0 c^2 {u_ {i+1/2}^{n+1/2} -U_ {i-1/2}^{n+1/2}} {\ delta x}} \ right) \ end {gather} The difference formula (3) is as follows. \ begin {gather} \ Frac {p_i^n-p_ {i-1}^n} {\ delta x} =-\ rho_0 \ frac {u_ {i-1/2}^{n+1/2} -U_ {i-1 /2}^{n-1/2}} {\ Delta t} \ nonumber \\ [12pt] u_ {i-1/2}^{n+1/2} = u_ {i-1/2}^{n-1/2}-\ flac {\ delta t} {\ rho_0} \ leaft (leaft ( \ Frac {p_i^n-p_ {i-1}^n} {\ delta x} \ right) \ end {gather} The symbols used in the above formula are as follows.
$ P $: Sound pressure
$ U $: $ x $ particle speed
$ ∆ t $: Ingredient width of time
$ ∆ x $: Ingredients in the space
$ I $: Lattice position in the space area
$ n $: Time steps
1/2 in the above formula means that the sound pressure and air particle speed in the spatial grid and time steps are halved.
If you explain this mathematical model in an easy -to -understand manner, it is a way to predict the sound pressure and particle speed of the next time by using the previous time and particle speed.
When the above formula is extended to two -dimensional, the sound pressure and particle speed are distributed as shown in FIG. The sound pressure of the next time is estimated from the sound pressure of the current time and the particle speed on the top, bottom, left and right, and the particulate speed of the next time is in the same direction as the current particle speed (left and right in the $ x $ direction, $ and $. In the Y $ direction, I guess from the upper and lower sound pressure).
Figure 5 Space arrangement of sound pressure and air particle speed
The formulas are as follows. \ begin {gather} P_ {i, j}^{n+1} = p_ {i, j}^n - \ delta t \ rho_0 c^2 \ leaft ( \ Frac {u_ {i+1/2, j}^{n+1/2} -U_ {i-1/2, j}^{n+1/2}} {\ delta x}+\ frac {frac { v_ {i, j+1/2}^{n+1/2} -v_V_ {i, j-1/2}^{n+1/2}} {\ delta y} \ right) \ end {Gather} \ begin {gather} u_ {i-1/2, j}^{n+1/2} = u_ {i-1/2, j}^{n-1/2}-\ frac {\ delta t} {\ rho_0} \ Left ( \ Frac {p_ {i, j}^n-p_ {i-1, j}^n} {\ delta x} \ right) \ end {Gather} \ begin {gather} v_ {i, j-1/2}^{n+1/2} = v_ {i, j-1/2}^{n-1/2}-\ frac {\ delta t} {\ rho_0} \ Left ( \ Frac {p_ {i, j}^n-p_ {i, j-1}^n} {\ delta y} \ right) \ end {Gather} The symbols used in the above formula are as follows.
$ P $: Sound pressure
$ U $, $ V $: $ x $ and $ Y $ particle speed
$ ∆ t $: Ingredient width of time
$ ∆ x $, $ ∆ y $: The chopped width of the space
$ I $, $ J $: Lattice position in space area
$ n $: Time steps
The above is the simplest differentiation method of the FDTD method. However, an error occurs by differentiating the differential equation. There are many calculation methods with higher accuracy (decrease in errors), but as an example, see the spatial wider range of sound pressure and air particle speed, and estimate the results of the next time (2). , 4) There is a method called law [3].
(2) Stable conditions
After formulating acoustic analysis using the FDTD method, set the space width and time width. In general, it is calculated over a long time to seek the results of acoustic analysis. There is an appropriate space width and time width set value so that the calculation by the FDTD method does not diverge, which is called a stable or CFL condition. In the case of one dimension, the stable conditions are shown in the following formula. \ begin {gather} C \ Delta T \ LEQ \ Delta x \ end {gather} In the case of two dimensions, the formula below.
\ begin {gather} c \ delta t \ leq \ frac {1} {\ sqrt { \ BIGL (\ Frac {1} {\ Delta x} \ bigr) ^2 + \ BIGL (\ Frac {1} {\ delta y} \ bigr) ^2 }} \ end {gather} To put it simply, the distance of information propagation by one time step must be shorter than the engraved width of the spatial lattice. In the case of two or three dimensions, it is necessary to shorten it. The FDTD method uses information in the current space and estimates the future a little ahead. If you set the time width large, you will probably guess the future far away, and there will be a gap with the reality. On the other hand, the space width setting is relatively free, but the higher the frequency you want to analyze, the smaller the space width must be set.
(3) Boundary conditions
When performing calculations, it is realistically impossible to assume an infinite space, so it is necessary to set the boundary.
There are two types of boundary conditions commonly used for acoustic analysis by the FDTD method.
(3-1) Rigid wall conditions
It is a condition that the sound waves are completely reflected at the boundary. It is often used in design of general rooms, reverberation rooms, and diffusion materials. It can be set by fixing the particle speed of the air on the boundary to 0.
(3-2) Complete absorption conditions
It is a condition that sound waves are completely absorbed at the boundary. This condition is used if you want to perform analysis in virtual and free sound fields. Completely absorption conditions are more important than rigid wall conditions, but the settings are also complicated. By using a complete absorption condition, even if the calculation area is reduced, the calculation result is not affected by the reflection of the boundary, and the memory and the calculation time can be saved. However, if the absorption condition is not set correctly, an unnecessary reflection will occur and the calculation result will occur. Here are three types of complete absorption conditions that are often used here:
(3-2-1) MUR primary absorption boundary conditions [4] \ begin {gather} P_ {0, j}^{n + 1} = p_ {1, j}^{n} + \ leaft ( \ Frac {C-\ Delta X} {C+\ Delta x} \ Right) \ Left ( P_ {1, j}^{n+1} -p_ {0, j}^{n} \ right) \ end {gather} \ begin {gather} P_ {i, 0}^{n + 1} = p_ {i, 1}^{n} + \ leaft ( \ Frac {C-\ Delta Y} {C+\ Delta y} \ Right) \ Left ( P_ {i, 1}^{n+1} -p_ {i, 0}^{n} \ right) \ end {gather} Equation 11 and 12 are examples of the boundary of $ (x = 0) $ (x = 0) $ and $ (y = 0) $. This is the simplest formula among the frequently used absorption boundary conditions. In equation 11, to calculate the sound pressure on the left border, only the sound pressure on the side of it is used. Similarly, to calculate the sound pressure at the lower boundary in Equation 12, use only one higher sound pressure. For this reason, the sound waves that enter vertically on the boundary are completely absorbed, but only the vertical incident incidental incident is gone, and the horizontal components remain as a reflector.
(3-2-2) HIGDON secondary absorption conditions [5]
If you want to reduce the reflector from the horizontal component, use the secondary absorption boundary conditions shown in the following formula.
\ begin {Equation} \ begin {Split} P_ {0, j}^{n+1} = 2 \ left (\ frax \ frac \ delta t- \ delta x} {c \ delta t+\ delta x \ right) \ left (p_ {1, j} ^{n+1} -p_ {0, j}^{n} \ right)+\ leaft (\ frac \ delta t- \ delta x} {c \ delta t+\ delta x} \ right)^2 \ left (p_ {2, j}^{n+1} -2p_ {1, j}^{n}+p_ {0, j}^{n-1} \ right)-\\ 2τ \ leaft (\ Frac {C \ Delta T- \ Delta X} {C \ Delta T+\ Delta X \ Right) \ Left (p_ {0, j}^{n} -p_ {1, j}^{{{{{{{{{{{{{{{ n-1} \ right)+2 {1, j}^{n-1}-τ^2p_ {0, j}^{n-1} \ END {Split} \ end {Equation} \ begin {Equation} \ begin {Split} P_ {i, 0}^{n+1} = 2 \ leaft (\ flac {c \ delta t- \ delta y} {c \ delta t+\ delta y} \ Right) \ leaft (p_ {i, 1} ^{n+1} -p_ {i, 0}^{n} \ right)+\ leaft (\ frax {c \ delta t- \ delta y} {c \ delta t+\ delta y} \ right)^2 \ Left (p_ {i, 2}^{n+1} -2p_ {i, 1}^{n}+p_ {i, 0}^{n-1} \ right)-\\ 2τ \ leaft (\ frac {c \ delta t- \ delta y} {c \ delta t+\ delta y} \ right) \ leaft (p_ {i, 0}^{n} -P_ {i, 1}^{{{{{{{{{{{{{{{{{{ n-1} \ right)+2 {p_ {i, 1}^{n-1}-τ^2p_ {i, 0}^{n-1} \ END {Split} \ end {Equation} Here is $ τ = 0.995 $.
The above two boundary conditions are ways to mathematically erase reflective sounds.
(3-2-3) Complete conformity layer (Perfect Matched Layer: PML) [6]
It is a condition for further suppressing the reflected sound.
This boundary condition is a concept similar to the wedge -shaped sound absorbing material in the 響 響 room. The reflection sound is erased by completely matching the incident sound and the acoustic impedance on the border. However, like the sound absorbing layer of the 響 室 room, the sound absorbing effect under PML conditions depends on the number of border layers. You can get a good sound absorbing effect if you set the border layer thickly, but you need to increase the memory and calculation time accordingly. As an example, there is a report that the establishment of a 20 -layer sound absorbing layer can get a sufficient sound absorption effect. [3]
(4) Sound source
There is no particular restriction on setting the sound source using the FDTD method. However, if you set a sound source that does not exist in reality, you will get meaningless results. Sign wave and Gaussian Pulse sound source are used as sound sources.
Sign wave is used to observe a single frequency. By observing the calculation process of the frequency, you can check whether the stable and boundary conditions are appropriate.
Gaussian Pulse is a signal containing a wide frequency band. Since the sound pressure distribution of the space is only required as the initial condition, the continuous calculation in the time area is not required. The setting method of the sound source is as follows.
\ begin {gather} P (x, y) = e^{-\ left (\ flac {r} {r} \ right)^2} \ end {gather} \ begin {gather} r = \ sqrt {\ left (x-x_s \ right)^2+\ left (y-y_s \ right)^2} {, {,} {, {n \ delta h} \ nonumber \ end {gather} The symbols used in the above formula are as follows.
$ P $: Sound pressure
$ e $: Napier number
$ x $, $ Y $: Position in space area
$ x_s $, $ y_s $: The center of the sound source in the space area
$ N $: Control the arbitrary constant and frequency bandwidth contained in pulse
$ ∆h $: Space width (when space lattice is a square)
There are many other sound sources setting methods, but even if the equation of the sound source is continuous, it is regarded as discontinued (such as a narrow pulse) depending on the space and time width setting, and an error occurs in the analysis. Be careful.
Above, we have introduced four basic elements in the FDTD method. With some programming skills, simple sound analysis can be realized. As mentioned above, the calculation result in the spatial area can be used as a video, and the calculation result in the time area is output as a WAV file, which can be heard. It is a really useful numerical analysis method.
Research is still ongoing for each element. Its purpose is high precision and error reduction. With the progress of the calculator, speeding up is one of the research themes.
The audio analysis by this FDTD method is a major applied field, considering the spread sound field in the indoor space and the prediction of sound transmission in the open space. By expanding the mathematical model of the FDTD method, it is possible to analyze the porous material and the sound wave propagation in solids, and it can also analyze it with air. In other words, it is also possible to predict the sound analysis when there is a sound absorbing material in the room and the sound insulation effect between the next room by wall material.
The disadvantage of acoustic analysis by the FDTD method is that the calculation time and the memory used are enormous, so it is not yet a widely used acoustic analysis method. However, in recent years, high -speed calculations have been realized in many cases. Therefore, I hope that the acoustic analysis by the FDTD method will be available in more situations in the future.
References
[1] Ayumi Ushiyama, Hiroshi Nagatomo, Shinichi Sakamoto, Hideki Tachibana, "Sound field forecast of Komaba Convention Hall by 3D FDTD analysis," Japanese Society of Society Lecture Papers, Pp.983-984, 2005.
[2] Makoto Yokota, Shinichi Sakamoto, Hideki Tachibana, Seiko Ishii, "Numerical analysis and hearing of the squealing dragon phenomenon by the FDTD method," Japanese Society of Architecture Society Environmental Paper Collection, Vol. 73, No. 629, pp. 849-856, 2008.
https://doi.org/10.3130/aije.73.849
[3] TETSUYA SAKUMA, SHINICHI SAKAMOTO, TORU OTSURU, Computational Simulation in Architectual and Environmental Acoustics, Springer, 2014.
https://doi.org/10.1007/978-4-431-54454-8
[4] G. MUR, "Absorbing Boundary Conditions for the Finite-Difference Approximation of the time-domain Electromagnetic-Field Equations," in IEEE Transactions on Electromagnetic Compatibility, Vol. EMC-23, NO. 4, pp. 377-382, 1981.
https://doi.org/10.1109/TEMC.1981.303970
[5] Robert L. HIGDON, "Absorbing Boundary Conditions for Difference Approximations to the Multi-Dimensional Wave Equation," Mathematics of Computation, Vol. 47, No. 176, pp. 437-459, 1986.
https://doi.org/10.2307/2008166
[6] J.-P. BERENGER, “A Perfectly Matched Layer for the Absorption of Electromagnetic Waves,” Journal of Computational Physics, Vol. 114, No. 2, pp. 185–200, 1994.
https://doi.org/10.1006/jcph.1994.1159
Seki
First, I will briefly introduce how the FEM method and the BEM method are applied to speaker development. Since the FEM method can efficiently analyze solids and fluids, it can verify how the sound generated by the combination of the shape of the diaphragm and the material changes (see Fig. 1 (A). ) Since the BEM method can calculate the acoustic characteristics on the boundary, the orientation frequency characteristics of the speaker can be verified by calculating the acoustic characteristics of 1m or 2m, mainly speakers (Fig. 1 (B)). )reference).
The FDTD method is still rarely used for speaker development, and is widely used in the architectural sound field that emphasizes sound wave propagation [1, 2].
The FDTD method is a method that is originally used to analyze electromagnetic waves. Based on the Maxwell equation, the time change of electromagnetic distribution in the observation space can be calculated using a relevant formula between magnetic fields and electric fields. The same approach allows the sound analysis (Note 2) by the FDTD method for the same wave phenomenon.
Note 1: Standing state does not change depending on time.
Note 2: Acouter analysis calculates the sound wave shape at arbitrary points in the space, and also requires each sound indicator (frequency characteristics, reverberation time, sound insulation performance, etc.).
The interesting thing about the sound field analysis by the FDTD method is that it can seek phenomena such as reflection, diffraction, and interference when the sound waves arrive in an obstacle. By visualizing this time change as a video, you can clearly indicate the state of sound wave propagation.
Here is an example of sound wave propagation analysis. Figure 2 and 3 assume a square meter spatial area, and set the area around it, such as a free sound field, is set. In Figure 2, when an impulse sound source is given to the center of the space, it is possible to confirm that the sound wave is propagated as a circle centered on the sound source. Further, in Fig. 3, a square obstacle is placed in the lower right to observe the phenomenon of reflection and diffraction by obstacles.
Next, I will introduce the basic elements required for acoustic analysis using the FDTD method. It is the setting of the formula, stable conditions, boundary conditions and sound sources. There are more complicated methods for these elements, but all of them are applied to calculus and calculator theory, so I will omit it in this article.
(1) Formatization of the mathematical model of sound wave propagation and the FDTD method First, I will explain from the mechanism of sound wave propagation. When vocal cords and vibration plates of speakers begin to vibrate as the source of sound, the surrounding air particles begin to move due to this vibration. With this movement, the air density changes and pressure (sound pressure) changes (state 1). And due to pressure changes, the nearby air particles also begin to move (state 2). The sound waves are formed by the repetitive periodic movement of this air particles, and the sound waves are propagated from the sound source.
Express each of the above conditions with a simple one -dimensional mathematical model.
Condition 1: As shown in Fig. 4. A, the air particles are densely $ $ $ and speed $ U $ on the air grid to be observed, and the speed of $ U+du $ leaks at this instant lattice. The mass change rate can be expressed as follows. \ begin {gather} \ Frac {d \ rho^\ prime} {dt} dx = \ rho_0 u -\ rho_0 (u + du) = - \ rho_0 du \ nonumber \\ [10pt] \ Frac {d \ rho^\ prime} {dt} = - \ rho_0 {du} {dx} \ end {gather} Here $ ρ '$ is the air density in the lattice. In addition, since there is a relationship between $ p = ρ'c^2 $ from the air condition equation, the equation (1) is (2). \ begin {gather} \ Frac {dp} {dt} = - \ rho_0 c^2 {du} {dx} \ end {gather} Here $ C $ is the speed of sound in the air.
Condition 2: If a pressure difference $ DP $ occurs on both sides of the air grid as shown in Fig. 4.B, the air particles in the lattice move at $ U $ in the second law of the exercise. The relationship between the pressure and speed is as follows. \ begin {gather} \ Frac {dp} {dx} = -\ rho_0 {du} {dt} {dt} \ end {gather}
Equation (2) and (3) are called continuous equations of sound wave propagation and exercise equation. The FDTD method can convert these two differential equations to differential equations and calculate sound wave propagation only in four arithmetic operations. The difference formula (2) is as follows. \ begin {gather} \ Frac {p_i^{n+1} -p_i^n} {\ delta t} = - \ rho_0 c^2 {u_ {i+1/2}^{n+1/2} -U_ {i {i {i_ {I_U_ {{{{{{{{{{{{{} -1/2}^{n+1/2}} {\ delta x} \ nonumber \\ [12pt] P_i^{n+1} = p_i^n - \ delta t \ Left ( \ rho_0 c^2 {u_ {i+1/2}^{n+1/2} -U_ {i-1/2}^{n+1/2}} {\ delta x}} \ right) \ end {gather} The difference formula (3) is as follows. \ begin {gather} \ Frac {p_i^n-p_ {i-1}^n} {\ delta x} =-\ rho_0 \ frac {u_ {i-1/2}^{n+1/2} -U_ {i-1 /2}^{n-1/2}} {\ Delta t} \ nonumber \\ [12pt] u_ {i-1/2}^{n+1/2} = u_ {i-1/2}^{n-1/2}-\ flac {\ delta t} {\ rho_0} \ leaft (leaft ( \ Frac {p_i^n-p_ {i-1}^n} {\ delta x} \ right) \ end {gather} The symbols used in the above formula are as follows.
$ P $: Sound pressure
$ U $: $ x $ particle speed
$ ∆ t $: Ingredient width of time
$ ∆ x $: Ingredients in the space
$ I $: Lattice position in the space area
$ n $: Time steps
1/2 in the above formula means that the sound pressure and air particle speed in the spatial grid and time steps are halved.
If you explain this mathematical model in an easy -to -understand manner, it is a way to predict the sound pressure and particle speed of the next time by using the previous time and particle speed.
When the above formula is extended to two -dimensional, the sound pressure and particle speed are distributed as shown in FIG. The sound pressure of the next time is estimated from the sound pressure of the current time and the particle speed on the top, bottom, left and right, and the particulate speed of the next time is in the same direction as the current particle speed (left and right in the $ x $ direction, $ and $. In the Y $ direction, I guess from the upper and lower sound pressure).
The formulas are as follows. \ begin {gather} P_ {i, j}^{n+1} = p_ {i, j}^n - \ delta t \ rho_0 c^2 \ leaft ( \ Frac {u_ {i+1/2, j}^{n+1/2} -U_ {i-1/2, j}^{n+1/2}} {\ delta x}+\ frac {frac { v_ {i, j+1/2}^{n+1/2} -v_V_ {i, j-1/2}^{n+1/2}} {\ delta y} \ right) \ end {Gather} \ begin {gather} u_ {i-1/2, j}^{n+1/2} = u_ {i-1/2, j}^{n-1/2}-\ frac {\ delta t} {\ rho_0} \ Left ( \ Frac {p_ {i, j}^n-p_ {i-1, j}^n} {\ delta x} \ right) \ end {Gather} \ begin {gather} v_ {i, j-1/2}^{n+1/2} = v_ {i, j-1/2}^{n-1/2}-\ frac {\ delta t} {\ rho_0} \ Left ( \ Frac {p_ {i, j}^n-p_ {i, j-1}^n} {\ delta y} \ right) \ end {Gather} The symbols used in the above formula are as follows.
$ P $: Sound pressure
$ U $, $ V $: $ x $ and $ Y $ particle speed
$ ∆ t $: Ingredient width of time
$ ∆ x $, $ ∆ y $: The chopped width of the space
$ I $, $ J $: Lattice position in space area
$ n $: Time steps
The above is the simplest differentiation method of the FDTD method. However, an error occurs by differentiating the differential equation. There are many calculation methods with higher accuracy (decrease in errors), but as an example, see the spatial wider range of sound pressure and air particle speed, and estimate the results of the next time (2). , 4) There is a method called law [3].
(2) Stable conditions
After formulating acoustic analysis using the FDTD method, set the space width and time width. In general, it is calculated over a long time to seek the results of acoustic analysis. There is an appropriate space width and time width set value so that the calculation by the FDTD method does not diverge, which is called a stable or CFL condition. In the case of one dimension, the stable conditions are shown in the following formula. \ begin {gather} C \ Delta T \ LEQ \ Delta x \ end {gather} In the case of two dimensions, the formula below.
\ begin {gather} c \ delta t \ leq \ frac {1} {\ sqrt { \ BIGL (\ Frac {1} {\ Delta x} \ bigr) ^2 + \ BIGL (\ Frac {1} {\ delta y} \ bigr) ^2 }} \ end {gather} To put it simply, the distance of information propagation by one time step must be shorter than the engraved width of the spatial lattice. In the case of two or three dimensions, it is necessary to shorten it. The FDTD method uses information in the current space and estimates the future a little ahead. If you set the time width large, you will probably guess the future far away, and there will be a gap with the reality. On the other hand, the space width setting is relatively free, but the higher the frequency you want to analyze, the smaller the space width must be set.
(3) Boundary conditions
When performing calculations, it is realistically impossible to assume an infinite space, so it is necessary to set the boundary.
There are two types of boundary conditions commonly used for acoustic analysis by the FDTD method.
(3-1) Rigid wall conditions
It is a condition that the sound waves are completely reflected at the boundary. It is often used in design of general rooms, reverberation rooms, and diffusion materials. It can be set by fixing the particle speed of the air on the boundary to 0.
(3-2) Complete absorption conditions
It is a condition that sound waves are completely absorbed at the boundary. This condition is used if you want to perform analysis in virtual and free sound fields. Completely absorption conditions are more important than rigid wall conditions, but the settings are also complicated. By using a complete absorption condition, even if the calculation area is reduced, the calculation result is not affected by the reflection of the boundary, and the memory and the calculation time can be saved. However, if the absorption condition is not set correctly, an unnecessary reflection will occur and the calculation result will occur. Here are three types of complete absorption conditions that are often used here:
(3-2-1) MUR primary absorption boundary conditions [4] \ begin {gather} P_ {0, j}^{n + 1} = p_ {1, j}^{n} + \ leaft ( \ Frac {C-\ Delta X} {C+\ Delta x} \ Right) \ Left ( P_ {1, j}^{n+1} -p_ {0, j}^{n} \ right) \ end {gather} \ begin {gather} P_ {i, 0}^{n + 1} = p_ {i, 1}^{n} + \ leaft ( \ Frac {C-\ Delta Y} {C+\ Delta y} \ Right) \ Left ( P_ {i, 1}^{n+1} -p_ {i, 0}^{n} \ right) \ end {gather} Equation 11 and 12 are examples of the boundary of $ (x = 0) $ (x = 0) $ and $ (y = 0) $. This is the simplest formula among the frequently used absorption boundary conditions. In equation 11, to calculate the sound pressure on the left border, only the sound pressure on the side of it is used. Similarly, to calculate the sound pressure at the lower boundary in Equation 12, use only one higher sound pressure. For this reason, the sound waves that enter vertically on the boundary are completely absorbed, but only the vertical incident incidental incident is gone, and the horizontal components remain as a reflector.
(3-2-2) HIGDON secondary absorption conditions [5]
If you want to reduce the reflector from the horizontal component, use the secondary absorption boundary conditions shown in the following formula.
\ begin {Equation} \ begin {Split} P_ {0, j}^{n+1} = 2 \ left (\ frax \ frac \ delta t- \ delta x} {c \ delta t+\ delta x \ right) \ left (p_ {1, j} ^{n+1} -p_ {0, j}^{n} \ right)+\ leaft (\ frac \ delta t- \ delta x} {c \ delta t+\ delta x} \ right)^2 \ left (p_ {2, j}^{n+1} -2p_ {1, j}^{n}+p_ {0, j}^{n-1} \ right)-\\ 2τ \ leaft (\ Frac {C \ Delta T- \ Delta X} {C \ Delta T+\ Delta X \ Right) \ Left (p_ {0, j}^{n} -p_ {1, j}^{{{{{{{{{{{{{{{ n-1} \ right)+2 {1, j}^{n-1}-τ^2p_ {0, j}^{n-1} \ END {Split} \ end {Equation} \ begin {Equation} \ begin {Split} P_ {i, 0}^{n+1} = 2 \ leaft (\ flac {c \ delta t- \ delta y} {c \ delta t+\ delta y} \ Right) \ leaft (p_ {i, 1} ^{n+1} -p_ {i, 0}^{n} \ right)+\ leaft (\ frax {c \ delta t- \ delta y} {c \ delta t+\ delta y} \ right)^2 \ Left (p_ {i, 2}^{n+1} -2p_ {i, 1}^{n}+p_ {i, 0}^{n-1} \ right)-\\ 2τ \ leaft (\ frac {c \ delta t- \ delta y} {c \ delta t+\ delta y} \ right) \ leaft (p_ {i, 0}^{n} -P_ {i, 1}^{{{{{{{{{{{{{{{{{{ n-1} \ right)+2 {p_ {i, 1}^{n-1}-τ^2p_ {i, 0}^{n-1} \ END {Split} \ end {Equation} Here is $ τ = 0.995 $.
The above two boundary conditions are ways to mathematically erase reflective sounds.
(3-2-3) Complete conformity layer (Perfect Matched Layer: PML) [6]
It is a condition for further suppressing the reflected sound.
This boundary condition is a concept similar to the wedge -shaped sound absorbing material in the 響 響 room. The reflection sound is erased by completely matching the incident sound and the acoustic impedance on the border. However, like the sound absorbing layer of the 響 室 room, the sound absorbing effect under PML conditions depends on the number of border layers. You can get a good sound absorbing effect if you set the border layer thickly, but you need to increase the memory and calculation time accordingly. As an example, there is a report that the establishment of a 20 -layer sound absorbing layer can get a sufficient sound absorption effect. [3]
(4) Sound source
There is no particular restriction on setting the sound source using the FDTD method. However, if you set a sound source that does not exist in reality, you will get meaningless results. Sign wave and Gaussian Pulse sound source are used as sound sources.
Sign wave is used to observe a single frequency. By observing the calculation process of the frequency, you can check whether the stable and boundary conditions are appropriate.
Gaussian Pulse is a signal containing a wide frequency band. Since the sound pressure distribution of the space is only required as the initial condition, the continuous calculation in the time area is not required. The setting method of the sound source is as follows.
\ begin {gather} P (x, y) = e^{-\ left (\ flac {r} {r} \ right)^2} \ end {gather} \ begin {gather} r = \ sqrt {\ left (x-x_s \ right)^2+\ left (y-y_s \ right)^2} {, {,} {, {n \ delta h} \ nonumber \ end {gather} The symbols used in the above formula are as follows.
$ P $: Sound pressure
$ e $: Napier number
$ x $, $ Y $: Position in space area
$ x_s $, $ y_s $: The center of the sound source in the space area
$ N $: Control the arbitrary constant and frequency bandwidth contained in pulse
$ ∆h $: Space width (when space lattice is a square)
There are many other sound sources setting methods, but even if the equation of the sound source is continuous, it is regarded as discontinued (such as a narrow pulse) depending on the space and time width setting, and an error occurs in the analysis. Be careful.
Above, we have introduced four basic elements in the FDTD method. With some programming skills, simple sound analysis can be realized. As mentioned above, the calculation result in the spatial area can be used as a video, and the calculation result in the time area is output as a WAV file, which can be heard. It is a really useful numerical analysis method.
Research is still ongoing for each element. Its purpose is high precision and error reduction. With the progress of the calculator, speeding up is one of the research themes.
The audio analysis by this FDTD method is a major applied field, considering the spread sound field in the indoor space and the prediction of sound transmission in the open space. By expanding the mathematical model of the FDTD method, it is possible to analyze the porous material and the sound wave propagation in solids, and it can also analyze it with air. In other words, it is also possible to predict the sound analysis when there is a sound absorbing material in the room and the sound insulation effect between the next room by wall material.
The disadvantage of acoustic analysis by the FDTD method is that the calculation time and the memory used are enormous, so it is not yet a widely used acoustic analysis method. However, in recent years, high -speed calculations have been realized in many cases. Therefore, I hope that the acoustic analysis by the FDTD method will be available in more situations in the future.
References
[1] Ayumi Ushiyama, Hiroshi Nagatomo, Shinichi Sakamoto, Hideki Tachibana, "Sound field forecast of Komaba Convention Hall by 3D FDTD analysis," Japanese Society of Society Lecture Papers, Pp.983-984, 2005.
[2] Makoto Yokota, Shinichi Sakamoto, Hideki Tachibana, Seiko Ishii, "Numerical analysis and hearing of the squealing dragon phenomenon by the FDTD method," Japanese Society of Architecture Society Environmental Paper Collection, Vol. 73, No. 629, pp. 849-856, 2008.
https://doi.org/10.3130/aije.73.849
[3] TETSUYA SAKUMA, SHINICHI SAKAMOTO, TORU OTSURU, Computational Simulation in Architectual and Environmental Acoustics, Springer, 2014.
https://doi.org/10.1007/978-4-431-54454-8
[4] G. MUR, "Absorbing Boundary Conditions for the Finite-Difference Approximation of the time-domain Electromagnetic-Field Equations," in IEEE Transactions on Electromagnetic Compatibility, Vol. EMC-23, NO. 4, pp. 377-382, 1981.
https://doi.org/10.1109/TEMC.1981.303970
[5] Robert L. HIGDON, "Absorbing Boundary Conditions for Difference Approximations to the Multi-Dimensional Wave Equation," Mathematics of Computation, Vol. 47, No. 176, pp. 437-459, 1986.
https://doi.org/10.2307/2008166
[6] J.-P. BERENGER, “A Perfectly Matched Layer for the Absorption of Electromagnetic Waves,” Journal of Computational Physics, Vol. 114, No. 2, pp. 185–200, 1994.
https://doi.org/10.1006/jcph.1994.1159
Seki
