CC BY
30Electronic Journal «Technical Acoustics» http://www .ejta.org
2008, 18
E. A. Navarro1, J. Segura2, R. Sanchis, A. Soriano
Institut de Robotica - University of Valencia, Poligon de la Coma s/n 46980-Paterna (Valencia), 1 enrique.navarro@uv.es, 2jaume.segura@uv.es
Solving 2D acoustic ducts and membranes by using FDTD method
Received 28.08.2008, published 27.11.2008
The paper describes an application of the FDTD (Finite Difference Time Domain) method, combined with the discrete Fourier transform (DFT) for solving Eigenvalue problems in acoustical waveguides and membranes of arbitrary cross-section. The governing acoustic equations are discretised in a two-dimensional domain, in which boundary conditions are defined over pressure or velocity. At cut-off, the cross-section of a waveguide acts as a two dimensional resonator. Therefore, the spectral response gives the cut-off frequencies of the waveguide or the resonant frequencies of the membrane. Once each frequency is known, the application of the discrete Fourier transform also provides the spatial distribution of the pressure-velocity modes in the cross-section of the waveguide or membrane.
Keywords: finite differences in time domain, Euler equations, waveguide, membrane.
INTRODUCTION
The numerical modeling of acoustical structures is a subject of interest for sound synthesis due to the characteristics of sound dynamics under different conditions. An interesting application field is the simulation of musical instruments which is of interest in order to evaluate the sound interaction between the exciting elements and the resonators accurately. When a sharp and short sound pulse excites, propagation in a certain structure gives rise to resonance and the propagation of sound can be decomposed into an infinite sum of "acoustical modes" [1]. In this paper we demonstrate the applicability of the Finite Difference Time Domain method (FDTD) to the analysis of waveguide modes and two-dimensional elastic media.
The use of acoustic waveguides in signal processing has widespread. As the boundary conditions are complex, it is difficult to use a simple model and various numerical methods are available in the literature to find cut-off frequencies and field mode distributions (finite element method, integral operator methods etc.) [2, 3, 4, 5]. These numerical methods often require large computer resources and a good deal of programming effort. The finite element method (FEM), which is the preferred one because of its accuracy and versatility, even gives spurious solutions that more recent formulations try to avoid [6] when applied to vector Eigenvalue problem.
Yee proposed the Finite Difference Time Domain Method (FDTD) [7], to study electromagnetic wave propagation. Taflove, Holland, and Umashankar, [8, 9, 10] extended the algorithm and the FDTD method was successfully applied to scattering and circuit problems in electromagnetism. Thus far, FDTD is the most versatile and attractive technique used in electromagnetic wave propagation. A selection of the FDTD literature can be found in reference [11].
The FDTD method discretizes the Maxwell curl's equations in space and time [7], establishing a set of difference equations for the six electromagnetic field components. These difference equations are second order accurate in space and time, in which the electric field components are obtained from known past magnetic field components at each time instant. Also, the magnetic field components are obtained from the electric components following a leapfrog technique. A similar technique was used for elastic wave propagation in which the electric and magnetic fields were substituted by stress and velocity [12, 13]. The application of FDTD to acoustic propagation can be considered a special case of the former technique for analyzing stress-velocity [14, 15] in elastic materials. However, in acoustics, FDTD stands for the leapfrog scheme used in the modeling of pressure and velocity.
The FDTD method was used in acoustic scattering [16, 17], and also to study room acoustics [18, 19]. A limitation of FDTD is the high computational cost and the application of boundary conditions. This cost is determined by the required spatial discretization, the cell size lower than XI10 to avoid dispersion and the curved boundaries that are modeled in a staircase way. To reduce the computational cost, Poorter and Botteldooren proposed a subgriding algorithm [20], trying to overcome the staircase boundaries by using quasi-Cartesian grids [21]. The simplicity of the FDTD equations offers a great flexibility, allowing the inclusion of viscosity and nonlinear effect [22, 23]. Finally, FDTD was applied to model the propagation of ultrasounds in a homogeneous thermoviscous fluid by combining the nonlinear absorbing wave equation with the thermal diffusion equation [24].
In the present work, the FDTD method is proposed to analyze the cut-off frequencies and the pressure-velocity distribution of acoustical modes in hollow acoustic waveguides and membranes. A numerical procedure that extends the FDTD method to the study of modes in waveguides and membranes of arbitrary cross-section is presented. This numerical formulation can be used to compute cut-off frequencies and velocity as well as and pressure-field distributions; however, a simple programming and fast computation are used. Also, the numerical formulation can be easily extended to analyze waveguides that have arbitrary shape and surface impedance, provided a suitable discrete form of the boundary conditions can be established [16]. On the other hand, the FDTD algorithm has several advantages over other finite difference schemes. First, the FDTD method uses centered spatial and temporal differences to approximate the derivatives. This provides second-order accuracy in both time and space derivatives compared to other first-order schemes. Second, there is no need for special treatment of the edges. Finally, the leapfrog algorithm does not represent dissipation and only has a small amount of dispersion, which is negligible, compared with the physical dispersion of the structure under study [25].
We present the numerical formulation and results obtained with the technique when applied to structures with known analytical solution: rectangular and circular shape. Numerical calculations are performed using different mesh dimensions to examine the sensitivity of the model. The procedure presented is easily extended to an elliptical pipe, in which a rigorous mathematical analysis would imply the use of the special Mathieu functions. Finally, a guitar-shaped waveguide is also analyzed and the results are compared to those obtained by an experimental procedure [1]. Examples of field patterns obtained for some of the first modes are also illustrated.
1. FORMULATION 1.1. The hollow rigid pipe
Consider a hollow rigid pipe with arbitrary cross-section. The inner material is considered to be air or any gas, being homogeneous and isotropic with a density p and sound propagation velocity c .
The governing acoustic equations in three dimensions in the time domain are
dp ~dt
2 ( dv dvy dv pc2
y- + ■
dx dy dz
V s J
dvx _ 1 dp
dt p dx
dVy _ 1 dp
dt p dy
dvz _ 1 dp
(1) (2) (3)
. (4)
dt p dz
wherep is the acoustic pressure and v = vxux + vyuy + vzuz is the displacement velocity.
Assuming translational symmetry of the pipe along the z-axis, we look for harmonic propagating waves along the z-axis. These are the named acoustical propagating modes along the pipe, these fields will have dependence in the z coordinate and time in the form:
H(x, y, z, t) = h(x, y)-cos{a>t ± f5z), (5)
where f is the propagation constant; co = 2f , f being the frequency of the harmonic wave, and H is any field component of above (1)-(4) equations. The z derivatives in equations (1)-(4) become
dH - ± jfH. (6)
At the cut-off frequencies (f = fc), the acoustic wave does not propagate. Therefore, we have propagation constant f = 0 and the equations (1)-(4) become:
dp Ht
2 ( dv dvy } pc2
y
--1--
^ dx d
dx dy
dvx _ 1 dp
dt p dx
dVy _ 1 dp
(7)
(8) (9)
dt p dy
c>v
= o. (10)
dt
To solve for the modes and their cut-off frequencies, these equations are defined over the cross section of the acoustical waveguide with the proper boundary conditions. If we have a rigid pipe, the boundary conditions on the rigid walls of the pipe are
V tP-n| s = 0, (11)
where n is the outward unit vector perpendicular to the rigid inner surface S of the pipe.
Equations (7), (8) and (9) with the boundary condition of equation (11) and harmonic time dependence lead to the Helmholtz equation, an elliptic equation that traditionally was solved by spectral methods, Finite Differences (FD) or Finite Element (FE) techniques.
Next, we describe an alternative approach following a time domain procedure: Let us assume that a time domain excitation is introduced in the cross-section of the waveguide. The equations (7), (8), (9) with the boundary condition (11) will provide the time domain evolution of three-field components, p, vx, and vy, on the transversal section of the pipe. This is a two-dimensional time domain problem that is only valid at the cut-off frequencies of the propagating modes in the pipe. Thus, after introducing a pulsed excitation, only harmonic time dependence signals are allowed inside the two-dimensional cross-section. These harmonic signals will have frequencies equal to the cutoff frequencies of the propagating modes. The spectral amplitude of the incident pulse provides the amplitude of these resonant waves. Therefore, the problem of finding cut-off frequencies and propagating modes in the three-dimensional pipe becomes two-dimensional and corresponds to finding the resonant frequencies and field distributions of velocity and pressure in the cross-section of the pipe, considered as a two-dimensional resonator.
After a time pulsed excitation is established in the cross-section of the pipe, the total field pressure can be expressed as a superposition of modes, whose cutoff frequencies are included in the spectrum of the time domain pulse introduced:
P( ^ ^ t) = Z cm '<t>.m (^ y)'exp(j 2^fmt) , (12)
m
where (¡)m (x, y) is the space distribution of the mode m in the cross-section of the pipe, fm is its cut-off frequency, and cm is a complex constant.
Equations (7)-(9) are discretized over the cross-section of the waveguide, establishing a mesh in which the field components are defined at each node in an alternate fashion. The unit cell used is shown in Figure 1.
Figure 1.
2D-FDTD cell used in the 2D- Eigenvalue acoustical problems
From the FDTD algorithm, the field values are computed in alternate time steps. The numerical program that solves the discretized equations (7)-(9) has a calculation loop, in which p and vx - vy, are calculated at alternative time steps in order to achieve a second order
accuracy in the time derivatives. Each p -vx-vy calculation is performed at a time increment
At, giving the time evolution of the field components p - vx-vy .
Particularly, the pressure at a given mesh point (i, j) gives a time series:
pn (i,j) = p(iAx, jAy, nAt), (13)
where n = 0, 1, 2,..., N-1, and N; At is the total simulation time.
The application of the discrete Fourier transform (DFT) to eq. (13) will result in a discrete function of frequency:
N -1
Fk (i,j) = 2 Pn (i, j)-exp(-7 2nfknAt),
(14)
n=0
where fk = k / N-At .
This expression can be rewritten by using eq. (12):
N-1
„ .. ^ , .. 1 - exp(-j2nfknAt) Fk (i, j) = 2cJn (i, j)-- k '
n=0 1 - eXP^.j2n(fm - fk ))
(15)
On relating discrete and continuous Fourier transforms [26], Fk (i, j) may be expressed as
1
F (i j) = v jc^n(i,j)^^_
k('1} 2nAt h L+(r / At))'
which is periodic with period N.
Thus, if the sampling theorem is verified [26], the final expression for Fk (i, j) is
(16)
Fk (i,j) = 2"
jCm0m (i,j)
2nAt(fm - fk)
where k = 0, 1, 2..., N/2, which gives the frequency spectrum at the sampling point (i, j)
(17)
The frequencies corresponding to the maximum of \Fk (i,j) | will approximately give the cut-off frequencies of the m modes (fm). The frequency resolution achievable by the DFT is
Af = 1/ N At, (18)
which is related to the number of time steps N and the cell dimensions Ax and Ay by the stability criterion [8]:
c
V
> At.
1 1 . (19)
- + -
Ax Ay
Apart from this factor, the accuracy to be gained in the spectrum calculation also depends on the space-time simulation process undertaken by the FDTD. The sampling intervals in both space and time satisfy boundary conditions, including the effects of stair-casing in the modeling of curved boundaries, and make an inevitable numerical round-off.
The field distribution for each mode can be calculated once the cut-off frequency is known. The mode shapes are obtained by performing the DFT of the desired field components at the frequency fm over the entire cross-section of the waveguide.
1.2. Membranes
The membranes are the simplest two-dimensional resonators. The membranes have fixed edges, have an area density a and a surface tension T, constant all over its surface in the plane x-y. When a small area dx-dy is displaced a small distance dz, the surface tension T acts to restore it to equilibrium and the equation of motion is a wave equation [1]
^ = c2 V2 z (20)
dt2
where c = yjT/a is the velocity for transverse wave propagation, and the operator V2 is defined over the two dimensional section of the membrane, i.e. in Cartesian coordinates V 2 =il + .d 2
dx2 dy 2
Equation (20) is mathematically identical to the wave equation for sound obtained by combining (7)-(10). Thus, the analysis of a membrane with fixed edges can be made by using the algorithm in section 2.1, including minor changes: z is formally substituted by the pressure p , and the velocity of sound in air is substituted by the velocity for transverse wave
propagation (c = ^T/a ). Additionally, there is a substantial change in the boundary conditions. Now at the edges of the membrane there is a Dirichlet homogeneous boundary condition over the pressure, which means that the edges are fixed (z = 0 at edges):
P L = (21)
Therefore, extension of the FDTD procedure to analyze vibrating membranes is straightforward. In the following section, the utility of the technique presented is demonstrated, in the modeling of a rectangular and a circular membrane, but the procedure can be used to simulate any arbitrary cross-section membrane.
2. APPLICATIONS
The method is tested by obtaining the cut-off spectrum of a hollow rectangular pipe and a circular pipe. In both cases, in order to excite a large number of modes a space distribution of unit impulse functions is used in the excitation: at t = 0, the field is zero everywhere except at some arbitrarily located points. Square cells with Ax = Ay = Ah are used and the time step is in the border of the stability criterion: At = 0.7-Ah/c .
To study the numerical sensitivity of the algorithm, a rectangular waveguide with a cross-sectional area of 6.0x4.0 m2 is analyzed with six different discretisations. The hollow rectangular pipe is modeled with 12x8, 24x16, 48x32, 96x64, 192x128 and 384x256 meshes using square cells and c = 340 m/s. The number of total time steps is increased in each mesh in order to obtain the same accuracy in the spectral response (18). A larger mesh has a smaller Ah and, thus a smaller At, then Af = 1/ N•At is kept always to 0.1 Hz by increasing N. The results obtained for the first five modes are shown in Table 1. Cut-off frequencies have insignificant variation with the meshes; this is because the cross-section of the waveguide is perfectly modeled in all six meshes. Therefore, the 12x8 mesh with 9714 time steps is accurate enough; it takes 55 s of CPU time on a Pentium II (350 MHz) personal computer to cover a spectrum up to 500 Hz.
Table 1. Convergence analysis in the rectangular cross-section waveguide
Cut-off frequencies obtained with FDTD for the first four pressure modes in a rectangular rigid waveguide with 6 mx4 m cross-section. Six different mesh schemes are used (Ax = Ay = Ah).
Deviation of the numerical results from the corresponding analytical data is also provided.
Frequencies are given in Hertz
f10 % f01 % f20 % f21 % f02 %
12x8 28.34 0.04 42.33 0.4 56.33 0.6 70.79 0.05 83.8 1.4
24x16 28.34 0.04 42.45 0.11 56.56 0.19 70.79 0.05 84.7 0.4
48x32 28.34 0.04 42.45 0.11 56.56 0.19 70.79 0.05 84.9 0.11
96x64 28.34 0.04 42.33 0.4 56.33 0.6 70.79 0.05 84.9 0.11
192x128 28.34 0.04 42.33 0.4 56.33 0.6 70.79 0.05 84.9 0.11
384x256 28.34 0.04 42.33 0.4 56.33 0.6 70.79 0.05 84.9 0.11
analytical 28.33 -- 42.50 -- 56.67 -- 70.83 -- 85.0 --
In order to increase the accuracy gained in the spectral response, a longer time simulation is needed (large N) for a given mesh, (18)—(19). However, to obtain any significant improvement in the achievable frequency resolution, a great number of time iterations are necessary, which appears to be undesirable with a large cost in CPU time.
An analysis of the numerical results has been done by considering a circular pipe, with radius r = 0.5 m, when the boundaries of the mesh do not agree exactly with the actual boundaries. The boundary is modeled with a staircase edge, since the mesh is built with square cells. Table 2 shows the calculated cut-off frequencies and their percentage deviation from the theoretical results for various discretisation schemes: r / Ah =5, 10, 20, 30, 40, 50, and 100, (Ax = Ay = Ah , and At = 0.7-Ah/c, c = 340 m/s). Again, the number of total time steps is increased with increasing mesh size in order to achieve the same accuracy in the spectral response. In the 20x20 mesh (r / Ah =10) the CPU time consumed for a spectrum up to 500 Hz was 552 s (N = 97142). In Table 2, results are found to converge with increasing grid resolution, which suggests that curved boundaries are not a limitation for the algorithm.
Table 2. Convergence analysis in the circular waveguide
Cut-off frequencies obtained with FDTD for pressure modes 11, 21, 01, and 31, in a circular rigid waveguide with a radius r = 0.5 m. Seven different mesh schemes are used (Ax = Ay = Ah ).
Deviation of the numerical results from the corresponding analytical data is also provided.
Frequencies are given in Hertz
r/Ah f11 % f21 % f01 % f31 %
5 184.6 7.3 294.4 10.9 403.1 2.8 497.4 9.4
10 191.8 3.6 314.0 5.0 402.4 3.0 429.6 5.5
20 193.7 2.7 317.7 3.9 406.7 1.9 435.0 4.3
30 195.4 1.8 322.6 2.4 409.3 1.3 441.1 2.9
40 196.6 1.2 324.6 1.8 411.0 0.9 445.5 2.0
50 196.7 1.2 324.9 1.7 411.1 0.9 446.0 1.9
100 198.1 0.5 328.1 0.7 413.1 0.4 450.6 0.9
Analytical 199.07 -- 330.53 -- 414.72 -- 454.53 --
The analysis of membranes is made by using the same programs, but changing the boundary condition. In waveguide analysis it was defined on the normal velocity to the waveguide wall (11). In membrane analysis the boundary condition is set over pressure at the
rim (21). The velocity for the transverse wave propagation (c = ^T/a ) is set equal to the
velocity of sound in previous simulation, c = 340 m/s. By making these changes, the physical model is modified and also the spectral response. Exactly the same meshes and discretizations in space and time are repeated for rectangular and circular membranes with the same N, to get the same sensitivity, 0.1 Hz, in the spectral response. Table 3 presents the numerical results for the resonant frequencies of the circular membrane compared to theoretical values. Results are found to converge with increasing grid resolution in the same way as in the analysis of the circular waveguide. The results for the rectangular membrane are not shown, its resonant frequencies are coincident with the cut-off frequencies of the rectangular cross-section pipe ( f„m) for n, m>0 [1].
Table 3. Convergence analysis in the circular membrane
Resonant frequencies obtained with FDTD for pressure modes 01, 12, 21, and 22, in a circular membrane with a radius r = 0.5 m. Seven different mesh schemes are used (Ax = Ay = Ah ). Deviation of the numerical results from the corresponding analytical data is also provided.
Frequencies are given in Hertz
r/Ah f01 % f12 % f21 % f22 %
5 240.2 7.7 381.0 8.1 514.9 7.3 623.4 4.3
10 250.0 3.8 398.4 4.0 536.0 3.5 572.9 4.1
20 255.0 2.1 406.0 2.1 544.2 2.1 584.5 2.2
30 257.1 1.3 409.6 1.2 549.1 1.2 589.9 1.3
40 258.0 0.9 411.0 0.9 551.3 0.8 592.0 0.9
50 258.3 0.8 411.5 0.7 551.9 0.7 592.7 0.8
100 259.3 0.4 413.1 0.4 553.8 0.3 595.1 0.4
analytical 260.3 -- 414.6 -- 555.7 -- 597.4 --
Finally, we calculate the modes of elliptical and guitar-shaped pipes. The analytical treatment of the elliptical shape is especially difficult because it involves the use of the Mathieu special functions. However, in the FDTD modeling we simply define the boundary in a loop of the FORTRAN code. Results for the spectral response are shown in Figure 2. The first modes of the elliptical waveguide are shown in Table 4 and plots for these are given in Figure 3.a-b-c-d.
Figure 2. Spectral response in the elliptical shaped pipe: a/Ah =40 (blue), a/Ah =200 (red)
Table 4. Convergence analysis in the elliptical waveguide
Cut-off frequencies obtained with FDTD for several pressure modes in an elliptical rigid waveguide with major axis a = 0.5 m, eccentricity e=0.75. Results are normalized to the resonant frequency of the even-11 mode. Three different mesh schemes are used with Ax = Ay = Ah . Deviation of the numerical results from the corresponding analytical data is also provided
a/Ah fo11/fe11 % fe21/fe11 % fo21/fe11 % fe31/fe11 % fo31/fe11 %
40 1.4439 1.2 1.7844 1.0 2.0464 1.4 2.5523 0.8 2.6276 3.6
60 1.4562 0.3 1.7965 0.3 2.0464 1.3 2.5523 0.8 2.6783 1.6
200 1.4586 0.2 1.7992 0.2 2.0639 0.3 2.5655 0.2 2.7080 0.5
analyt. 1.4610 -- 1.8024 -- 2.0695 -- 2.5717 -- 2.7222 --
(c) (d)
Figure 3. Pressure distribution for some modes of the elliptical shaped pipe: (a) even-11, (b)
odd-11, (c) even-31, and (d) odd-31
The modeling of the guitar-shaped waveguide provides information about the influence of the guitar shape in its spectral response, because transversal resonance is produced at lower frequencies than longitudinal resonance between top and bottom plates. Fig. 4 shows the pressure distribution of the first mode, and Table 5 shows the resonant frequencies obtained. The comparison is difficult, because it depends on the exact shape of the guitar. However, a larger discussion about the guitars will be the focus of investigation in the future. From the results obtained some conclusions are evident: The error in the frequency versus the cell size is linear in curved boundaries, in which a staircase approach is involved. By increasing the
mesh size (decreasing the cell size), the "numerical boundaries" converge to the actual boundaries. For each mesh we have a different numerical domain that approaches the actual domain. Therefore, the second order accuracy is lost in getting the resonant frequencies when staircase boundaries are involved [27]. However, when the "numerical boundaries" match the actual boundaries, second order accuracy is achieved in the field calculation. The error in the resonant frequencies can be associated to the resolution in the location of the maximum and nulls of the fields along the numerical domain, i.e. the rectangular waveguide (corresponding to Table 1). Results showing the second-order accuracy obtained are not shown but can be found in reference [27].
(a) (b) (c)
Figure 4. Pressure distribution of the mode A1 (a), A3 (b) and A4 (c)
Table 5. Cut-off frequencies of the first ten modes in a guitar cross-section pipe
In the second and third line of the table, the first six modes of the air cavity of the guitar Martin (models D28 and D35) from reference [1] can be seen
Mode 1 2 3 4 5 6 7 8 9 10
Frequency 126.8 371.1 493.4 543.5 581.3 667.6 736.2 749.6 778.9 825.1
Martin D28 121 (A0) 383 (A1) 504 (A2) 652 (A3) 722 (A4) 956 (A5)
Martin D35 118 (A0) 392 (A1) 512 (A2) 666 (A3) 730 (A4) 975 (A5)
CONCLUSIONS
An alternative, simple and efficient method for studying hollow acoustical waveguides and membranes has been described. The sensitivity of the model to the cell dimensions and staircase boundaries has been examined by comparing the numerical results derived for various grid sizes for rectangular and circular shapes. The staircase approximation for non-planar boundaries is found to influence the calculation of the resonant frequencies. The frequency error versus cell size becomes linear with a slope depending on the calculated mode. However, in Eigenvalue problems in which the numerical boundaries coincide with the actual boundaries, the second-order accuracy with the cell size is recovered.
ACKNOWLEDGEMENTS
This research has been supported by the Universitat de València.
REFERENCES
1. N. H. Fletcher and T. D. Rossing. The physics of Musical Instruments. SpringerVerlag, New York 1991.
2. M. Guerich and M. A. Hamdi. A numerical method for vibro-acoustic problems with incompatible finite element meshes using B-spline functions. J. Acoust. Soc. Am., vol. 105, no. 3, pp. 1682-1694, March 1999.
3. F. Fontana and D. Rocchesso. Physical Modeling of Membranes for Percussion Instruments. ACUSTICA-acta acustica, vol. 84, pp. 529-542, 1998.
4. J. Bretos, C. Santamaria and J. Alonso-Moral. Vibrational patterns and frequency responses of the free plates and box of a violin obtained by finite element analysis. J. Acoust. Soc. Am., vol. 105, no. 3, pp. 1942-1950, March 1999.
5. A-C. Hladky-Hennion and R. Bossut. Time analysis of immersed waveguides using the finite element method. J. Acoust. Soc. Am., vol. 104, no. 1, pp. 64-71, Jul. 1998.
6. T. Angkaew, M Masanori and N. Kumagai. Finite element analysis of waveguide modes: a novel approach that eliminates spurious modes. IEEE Trans. Microw. Theory Techniques, vol. MTT-35, (2), pp. 117-123, 1987.
7. K. S. Yee. Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media. IEEE Trans. Antenna Propagat., vol. AP-14, (5), pp. 302-307, 1966.
8. A. Taflove and M.E. Brodwin. Numerical solution of steady-state electromagnetic scattering problems using the time dependent Maxwell's equations. IEEE Trans. Microw. Theory Techniques, vol. MTT-23, (8), pp. 623-630, 1975.
9. R. Holland, L. Simpson and K. Kunz. Finite difference analysis of EMC coupling to lossy dielectric structures. IEEE Trans. Electromagnetic Compat., vol. EMC-22, pp. 203-209, 1980.
10. K. Umashankar and A. Taflove. A novel method to analyze electromagnetic scattering of complex objects. IEEE Trans. Electromagnetic Compat., vol. EMC-24, pp. 397-405, 1982.
11. Kurt L. Shlager and John B. Schneider. A selective survey of the Finite-Difference Time-Domain Literature. IEEE Antennas and Propagation Magazine, vol. 37, no.4, August 1995, pp.39-56.
12. R. Madariaga. Dynamics of an expanding circular fault. Bull. Seismol. Soc. Am., vol. 66, pp. 639-666, 1976.
13. J. Virieux. SH-wave propagation in heterogeneous media: Velocity-stress finite difference method. Geophysics, vol. 49, pp. 1933-1942, 1984.
14. Q-H. Liu, E. Schoen, F. Daube, C. Randall, H-L. Liu and P. Lee. A three-dimensional finite difference simulation of sonic logging. J. Acoust. Soc. Am., vol.100, no.1,
pp. 72-79, 1996.
15. Y-H. Chen, W. C. Chew and Q-H. Liu. A three-dimensional finite difference code for the modeling of sonic logging tools. J. Acoust. Soc. Am., vol. 103, no. 2, pp. 702712, 1998.
16. S. Wang. Finite-difference time-domain approach to underwater acoustic scattering problems. J. Acoust. Soc. Am., vol. 99, no. 4, pt. 1, pp. 1924-1931, Apr.1996.
17. F.D. Hastings, J.B. Schneider and S.L. Broschat. A finite-difference time-domain solution to scattering from a rough pressure-release surface. J. Acoust. Soc. Am., vol. 102, no. 6, pp. 3394-3400, 1997.
18. D. Botteldooren. Finite-difference time-domain simulation of low-frequency room acoustic problems. J. Acoust. Soc. Am., vol. 98, pp. 3302-3308, 1995.
19. J. LoVetri, D. Mardare and G. Soulodre. Modeling of the seat dip effect using the finite-difference time-domain method. J. Acoust. Soc. Am., vol. 100, pp. 2204-2212,
1996.
20. J. De Poorter and D. Botteldooren. Acoustical finite-difference time-domain simulations of subwavelength geometries. J. Acoust. Soc. Am., vol. 104, no. 3, pp. 1171-1177, 1998.
21. D. Botteldooren. Acoustical finite-difference time-domain simulation in quasi-Cartesian grids. J. Acoust. Soc. Am., vol. 95, no. 5, pp. 2313-2319, 1994.
22. D. Botteldooren. Vorticity and entropy boundary conditions for acoustical finite-difference time-domain simulations. J. Acoust. Soc. Am., vol. 102, pp. 170-178,
1997.
23. D. Botteldooren. Numerical model for moderately nonlinear sound propagation in three-dimensional structures. J. Acoust. Soc. Am., vol. 100, pp. 1357-1367, 1996.
24. I. M. Hallaj, R. O. Cleveland. FDTD simulation of finite-amplitude pressure and temperature fields for biomedical ultrasound. J. Acoust. Soc. Am., vol. 105, no. 5, pp. L7-L12, 1999.
25. X. Zhang and K. K. Mei. Time-domain finite difference approach to the calculation of the frequency-dependent characteristics of microstrip discontinuities. IEEE Trans., 1988, MTT-36, (12), pp. 1775-1787
26. E. O. Brigham. The fast Fourier transform and its applications. Prentice-Hall, 1988.
27. E. A. Navarro, N. T. Sangary and J. Litva. Some Considerations on the accuracy of the non-uniform FDTD method and its application to waveguide analysis when combined with the Perfect Matched Layer. IEEE Trans. Microwave Theory Tech., vol. 44, no. 7, pp. 1115-1124, July 1996.
