Finite difference time domain modelling of a point-excited elastic plate radiating into an acoustic cavity

Finite difference time domain (FDTD) models are developed to solve the vibroacoustic problem of a thin elastic plate undergoing point force excitation and radiating into an acoustic cavity. Vibroacoustic modelling using FDTD can be computationally expensive because structure-borne sound wavespeeds are relatively high and a fine spatial resolution is often required. In this paper a scaling approach is proposed and validated to overcome this problem through modifications to the geometry and physical properties. This allows much larger time steps to be used in the model which significantly reduces the computation time. Additional reductions in computation time are achieved by introducing an alternative approach to model the boundaries between the air and the solid media. Experimental validation is carried out using a thin metal plate inside a small reverberant room. The agreement between FDTD and measurements confirms the validity of both approaches as well as the FDTD implementation of a thin plate as a three-dimensional solid that can support multiple wave types. Below the lowest room mode, there are large spatial variations in the sound field within the cavity due to the radiating plate; this indicates the importance of having a validated FDTD model for low-frequency vibroacoustic problems.


I. INTRODUCTION
Prediction of low-frequency sound fields inside cavities with airborne sources is well-suited to finite difference time domain (FDTD) models (e.g., see Refs. 1 and 2).However, many noise control problems in engineering involve lowfrequency sound radiation from vibrating surfaces into acoustic cavities.Hence this paper concerns FDTD modelling of the interaction between the vibration field of a pointexcited thin elastic plate and the resulting sound field at low frequencies inside an acoustic cavity.Running such a model can be problematic because the general FDTD formulation 3 that uses a uniform FDTD grid requires a large number of calculation cells and small time steps.
In this paper a scaling approach is introduced for FDTD modelling in order to give significant computational advantages.Previous approaches to model fine geometric details embedded in large FDTD models have primarily been concerned with the use of non-uniform grids, parallelization of the FDTD computations, and the use of dimension-reduced models.Parallelization requires the splitting of routine calculations and memory over a number of central or graphics processing units to reduce computation times and has been implemented by a variety of authors in a number of research fields such as vibroacoustics, 3 electrodynamics, 4 and vibration. 5The use of non-uniform grids (a technique also referred to as sub-gridding) allocates a finer grid resolution to regions that require more detail whilst using coarser grid cells elsewhere. 4Subgridding has been used in acoustic FDTD problems 6,7 but does not seem to have been considered in vibroacoustics.
Dimension-reduced models use a two-dimensional grid to represent solid media rather than three-dimensional calculation cells.This results in significant memory savings and reduced computation times.Previous work used KirchhoffÀLove theory for a two-dimensional implementation of thin elastic plates 8 and MindlinÀReissner theory for two-dimensional implementation of thick elastic plates. 9hese approaches are more complex to implement than general viscoelastic explicit FDTD routines.To provide flexibility in modelling, this paper considers a general viscoelastic explicit formulation for a solid medium which allows all forms of wave motion.In practice, bending waves are usually of most interest for vibroacoustic problems due to their relatively efficient sound radiation.
This paper focuses on a point-excited elastic plate radiating into an acoustic cavity.Section II describes the theory used to create the FDTD model and introduces an alternative approach to modelling boundaries between acoustic and solid media.A scaling approach is also introduced to improve computational efficiency for explicit methods that are limited by the Courant condition (but avoids the need for sub-gridding schemes).This approach is valid for vibroacoustic problems where the structural response is dominated by bending wave motion.Section III describes the experimental setup used to validate the FDTD model.Section IV contains results and analysis from the comparison of measurements and FDTD.

II. THEORY A. FDTD model
Previous work by Toyoda and Takahashi 3 modelled a vibrating solid in an acoustic cavity using FDTD by a) Electronic mail: carl.hopkins@liv.ac.uk considering only velocity and stress fields.The viscoelastic field is described using a velocity field (vector) and a stress field (tensor).However, in this work a pressure field (scalar) is used for the acoustic medium.This makes it easier to integrate existing acoustic FDTD routines, particularly those implementing boundaries using a perfectly matched layer (PML). 10Velocity and stress fields are used to model the solid medium, whereas velocity and pressure fields are used to model wave propagation in the acoustic medium, including PML boundaries.Coupling between the acoustic medium and the solid medium is a sequential two-way process such that any change in the acoustic field is transferred to the vibration field of the solid medium and vice versa.This is implemented using the following procedure: first, updating pressure field nodes, second, updating stress field components of the solid medium, and third, updating all the velocity components for the air and solid media.This updating loop then repeats until the simulation has reached a prescribed value of time, upon which it is stopped.
Stress nodes in the solid medium that are located outside of the plate, and acoustic pressure nodes that are located inside the plate, are not considered for the calculations and need not be allocated memory.

Acoustic wave propagation
The set of acoustic update equations to be used in the FDTD routine are described by Kunz. 11 These are the equation of motion for air particles, and continuity, where p is the acoustic pressure, v i is the velocity vector in the ith coordinate corresponding to the x-, y-, or z-directions, q o is the density of air and c is the speed of sound in air.

Structure-borne sound propagation
A general formulation is used for the solid medium that allows all types of wave motion.Momentum equations are then formulated in terms of the relationship between the stress tensor acting on an element of the solid material and the corresponding velocity of that element. 3The basic constitutive model assumed for the solid medium is a linear viscoelastic Voigt-element. 12Time domain equations using Einstein summation notation that describe the wave propagation are given by 3,13 q @v i @t ¼ @r ji @x j À bv i ; (3) where q is the material density, r ij is the stress tensor, d ij is the Kronecker delta, k and l denote the Lam e constants (which can be derived from other material properties 13 ), b is a damping constant proportional to velocity, and v and c denote the viscosity damping constants.The variables e ij and e kk are defined by In order to implement the FDTD routine, Eqs. ( 3) and ( 4) are discretized to allow their explicit numerical solution.
As an example of a set of discretized equations, consider the simple case of a two-dimensional space lattice (see Fig. 1) and consider the stress and velocity components along the zdirection as described by Eqs. ( 7) and ( 8): The lattice in Fig. 1 shows how the two-dimensional components of stress (represented by r and s for the normal and tangential components, respectively) and velocity (whose nodes are represented by double arrows) are placed at different positions forming a staggered grid.However, the normal stress components, r xx and r zz , are placed at the same position.In addition to the spatial offset, there is also an offset in time, where the stress components are calculated at different times to the velocity components.
Based on the space lattice shown in Fig. 1, Eq. ( 7) is discretized as follows: Dt q h D i r n xz j iÀ0:5;jþ0:5 þD j r n zz j i;j À bv nÀ0:5 z j i;jþ0:5 i ; where D i denotes the forward difference for spatial point i and n denotes the time step.For example, In order to discretize Eq. ( 8), the mixed derivatives can be reversed using Clairaut's theorem: where a z represents the acceleration in the z direction.In discrete form it is calculated from v z using Using this variable substitution, Eq. ( 8) is discretized as r nþ1 zz j i;j ¼r n zz j i;j þDt kD i v nþ0:5 x j iÀ0:5;j h þðkþ2lÞD j v nþ0:5 z j i;jÀ0:5 þvD i a nþ0:5 x j iÀ0:5;j þðvþ2cÞD j a nþ0:5 z j i;jÀ0:5 i : Equations ( 9) and ( 14) essentially form the update equations for v z and r zz , respectively.The nature of the method is that the value of a variable at a given instant of time is calculated directly from the values obtained at previous instants of time.Note that in order to transform Eqs. ( 9) and ( 14) into a form that can be implemented using a programming language, the indexes of the field array variables must be integers and cannot be fractions.
For the vibroacoustic problem considered in this paper, the variables are arranged in three-dimensional space according to a lattice described by Schroeder et al. 14

Vibration source
The simulations use a "hard" vibration source as described by Schneider. 15For this source the normal stress component in the z-direction is assigned a time-dependent force, F(t), which is converted into r zz ðtÞ through division by the area of a single horizontal grid cell, DxDy.The vibration source node follows the driving function irrespective of the state of its neighbour nodes.

Damping
Models are intended for the low-frequency range; hence, damping is not considered for wave motion in air.However, it is essential to include damping for the solid medium.All dissipation of mechanical energy occurs with internal damping due to the viscoelastic nature of the plate material which is incorporated using the approach of Toyoda et al. 21Damping constants b, c, and v result in a frequency-dependent loss factor, similar to systems characterized by Rayleigh damping.The constant b results in damping, which is proportional to velocity.The frequency-dependence of the b loss factor curve is inversely proportional to frequency.Constants c and v give equivalent frequency characteristics; hence it is possible to consider only c, and set v to zero whilst still obtaining a general Rayleigh damping profile.Using the constants b and c, the frequency-dependent loss factor is then given by 22 where f is the frequency and g indicates the loss factor.The frequency-dependent loss factor of the plate is determined from measurements.The damping coefficients are calculated such that the loss factor used in FDTD follows the measured loss factor as closely as possible; this is described in Sec.III B.

Numerical grid positions
Sound pressure is calculated at discrete points along uniform horizontal and vertical grids.The geometrical positioning of these grids matches the grids used in the experimental validation.The pressure impulse response at the grid positions is Fourier transformed and divided by the spectrum of the driving function FðxÞ, in order to obtain a transfer function of pressure-to-force.

Boundary conditions for the acoustic medium
For the room, the acoustic boundary conditions are frequency-independent and are described in detail by Yokota et al. 16

Boundary conditions around the edges of a plate
For the edges of the plate, the implementation of simply supported boundaries aims to approximate the following conditions: 17 where w denotes displacement and M x and M y indicate bending moments along the x-and y-directions, respectively.
To implement simply supported boundaries using the general viscoelastic FDTD formulation, only the kinematic condition w ¼ 0 needs to be specified.This is approximately carried out by assigning a value of zero to the vertical velocities that are located on the mid-plane around the plate edges as shown in Fig. 2. In this diagram the velocity nodes of the air particles are represented by a single arrow whereas the velocity of the solid medium is represented by a double arrow.Note that the approach described here differs from the implementation used for dimension-reduced models, 8 which require both displacement and bending moment conditions.

Boundary conditions at the interface between the plate and the surrounding air
For velocity nodes at the boundary between air and solid media, Toyoda et al. 18 split the velocity update equation into two equations involving a forward difference and a backward difference.These equations are combined to form a new equation, where the space step across the boundary is divided by a factor of 2 and the density at the boundary between the two media is averaged.This type of boundary condition is referred to here as the "standard approach."To improve computational efficiency, the implementation in this paper considers the update equations for the velocity nodes that lie on the boundaries to have the same form as the other solid medium velocity update equation [Eq.( 9)] for which the density equals that of the actual solid.However, both the pressure and stress field are modelled; hence the velocity update equation for a boundary node must include both pressure and stress terms, as indicated in Eq. ( 17): It is assumed that compression corresponds to a positive pressure increment, whereas in terms of stress tensors the same compression is assumed to be a negative stress increment (r xx ¼ r yy ¼ r zz ¼ -p). 19This sign convention is used in Eq. ( 17).The shear stress nodes adjacent to the boundary are set to zero because air is assumed to be an inviscid medium.In this paper, the approach to implementing boundary conditions will be referred to as a "simplified approach" because the implementation requires fewer calculations than the standard approach.For convenience the derivation considers a plate lying in a Cartesian coordinate plane, although it is feasible (but more complex) to consider other plate orientations.
The standard and simplified approaches result in the pattern shown in Fig. 3.Note that the principle applies to all nodes of the solid medium that are adjacent to the surrounding medium, not just the upper and lower surfaces of the solid medium.In this lattice diagram, the simplified approach results in boundary conditions for the plate which appear to be "ragged," although in terms of its eigenfrequencies the plate behaves as if it has smooth edges.Note that the eigenfrequencies correspond to a different set of physical dimensions as might be expected from the number of nodes in the solid medium.The physical dimensions along a given direction using the simplified boundary conditions are referred to as "effective" dimensions, which are L x;eff ; L y;eff , and L z;eff .Considering the number of normal stress nodes assigned to the plate, or by considering the standard approach, the physical dimensions that are to be expected, are referred to as "expected" dimensions, L x; exp ; L y; exp ; and L z; exp .In Fig. 3, the spatial offset along the ith Cartesian direction between L i;eff and L i; exp is denoted by d i;boundary .L i;eff can be obtained from L i; exp via the relation L i;eff ¼ L i; exp À 2d i;boundary .Hence, it is necessary to quantify d i;boundary in order to predict L i;eff from L i; exp and to be able to design plates with a prescribed set of dimensions.
To quantify the value of d i;boundary along a given direction, i, of the plate, a number of numerical tests are now carried out.These tests are based on the assumption that the mismatch between the expected and effective dimensions along a given direction must vanish as the number of calculation cells in that direction is increased.Therefore, tests are carried out with the minimum possible number of nodes along the thickness direction of the plate (two normal stress nodes) and the largest possible number of nodes along the other two lateral directions.The difference between the eigenfrequencies obtained from the FDTD model using the simplified approach and the eigenfrequencies corresponding to a plate with the expected dimensions given by Eq. ( 19) is primarily due to the difference between the expected and the effective thickness of the plate.Using this numerical approach it is only possible to identify an approximate value for d i;boundary ; however, this approximation becomes more accurate as the number of lateral calculation cells increases.
To carry out these tests, consider a simply supported plate in the xy plane for which the analytical equation describing the eigenfrequencies f p;q for thin plate bending wave motion is given by 20 where L z is the plate thickness, c L is the quasi-longitudinal wavespeed, and L x and L y are the lengths of the plate along the x-and y-directions, respectively.The eigenfrequencies of a plate with expected dimensions are calculated as follows: and the eigenfrequencies of a plate with effective dimensions are given by As the number of lateral nodes along the x-and y-directions is increased relative to the z-direction, the accuracy is increased for the two approximations: L x;eff % L x; exp and L y;eff % L y; exp .With the approximations for the expected and effective lengths across the lateral directions of the plate, the expected and effective eigenfrequencies are then related to the expected and effective thicknesses of the thin plate by Since L z;eff ¼ L z; exp À 2d z;boundary it can be shown that To allow bending wave motion, the minimum number of normal stress nodes across the thickness (z-direction) of the plate is two.This allows the lower stress node to be strained and the upper node to be under compression, and vice versa.Two thickness nodes are assumed in the following analysis because there is no computational advantage in using higher numbers of thickness nodes, but the same process can be followed if this is required.Hence, L z; exp ¼ 2Dz and the following equation relates d z;boundary to the spatial resolution along the z-direction: The result in Eq. ( 23) can also be used to relate d i;boundary to the spatial resolution along the ith direction.
In principle, Eq. ( 23) applies to any plate mode; however, the numerical tests use the fundamental mode (p ¼ q ¼ 1) to determine the effective thickness because higher modes are affected by numerical errors such as spatial discretization and numerical dispersion.
The numerical tests are carried out using undamped plates (arbitrary material properties) with two normal stress nodes along the thickness direction (for example, see Fig. 3 where i corresponds to the z-direction) and a varying number of normal stress nodes along each of the horizontal x-and y-directions.The spatial resolution is set to 0.025 m in all directions so the expected plate thickness h exp is 0.05 m.The results are given in Table I in terms of the ratio f eff 1;1 =f exp 1;1 .Assuming that the mismatch between the expected and effective lengths along a given direction vanishes as the number of nodes is increased in that direction, the ratio of 0.70 obtained using 120 stress nodes can be considered to be the most accurate; hence the value of d z;boundary is calculated to be d z;boundary ¼ ð1 À 0:7ÞDz ¼ 0:3Dz.Therefore, the relation between d i;boundary and the spatial resolution along the ith direction, denoted by x i , is also given by The next step is to define a methodology using the simplified boundaries approach to model a generic plate with a prescribed set of dimensions.In order to model a plate with n normal stress nodes and a side length of L x;eff along the xdirection, it is necessary to set the spatial resolution, Dx, to L x;eff ¼ nDx À 2 Â 0:3Dx () Dx ¼ L x;eff =ðn À 0:6Þ.It can be seen that the spatial resolution which is required along that direction is defined by Dx ¼ L x =ðn À 0:6Þ rather than Dx ¼ L x =n, where n is the number of normal stress nodes along a given direction.If n is set to two, which represents the number of stress nodes along the thickness direction, the space step along the thickness direction needed to implement the simplified boundary conditions approach is around 40% larger than that required using the approach described by Toyoda et al. 18 This larger space step provides significant computational benefits because the FDTD time step will be larger; hence the number of iterations required to reach a given time interval will be reduced.

C. Scaling of vibroacoustic fields
The Courant condition determines the maximum possible value for the time step in the FDTD model according to the following inequality: 11 Dt C where C is the highest phase velocity of any wave motion within the frequency range of excitation and Dx; Dy and Dz are the spatial resolutions in the x-, y-, and z-directions, respectively.If Dt does not satisfy the inequality in Eq. ( 24), the solution becomes unstable, i.e., the response will be unbounded.In practice, other factors such as damping and boundary conditions can also give rise to unstable solutions and it is not always possible to identify whether waves that have the highest phase velocity have actually been excited.Hence the aim is to find the largest value of Dt that provides a stable solution over the time period of the FDTD simulation.
In contrast to room acoustic simulations, it can be computationally expensive to run a large vibroacoustic model with a fine spatial resolution because of the fact that structure-borne sound wavespeeds are relatively high.In such cases, it is proposed to use an alternative formulation for the vibroacoustic problem to yield much faster results than with standard FDTD. 3The main issue with vibroacoustic models is that if the spatial resolution is kept constant throughout the domain, a large number of cells are required to represent the whole domain.With a fine spatial resolution, the required time step will be very small as a consequence of the Courant condition.To overcome this problem, the proposal is to model a larger structure that has the same vibration characteristics as the actual structure and couple it to the acoustic medium because a coarse spatial resolution will result in a larger time step.Once the equivalent structure is identified and processed, the results can be scaled back to represent the actual structure.Assuming the thickness direction of the plate is coincident with the z-direction (vertical direction), the following steps are used to scale the vibroacoustic model: (1) A scaling factor s > 1 is chosen and a plate with the same eigenfrequencies as the actual plate is identified where the side dimensions of the scaled plate are L x 0 ¼ sL x and L y 0 ¼ sL y ; respectively.In order to obtain the same bending eigenfrequencies [given by Eq. ( 18)], the thickness of the scaled plate is The spatial resolution of the scaled problem is then dictated by the dimensions of the scaled plate to give Dx 0 ¼ sDx, Dy 0 ¼ sDy; and Dz 0 ¼ s 2 Dz.In addition to scaling the plate, the x-, y-, and z-dimensions of the cavity also need to be scaled up by a factor of s to match the scaled x-and y-dimensions of the plate and to maintain the eigenfrequencies of the cavity.The scaled acoustic cavity with dimensions L x 0 ¼ sL x ; L y 0 ¼ sL y ; and L z 0 ¼ sL z is modelled using the aforementioned spatial resolutions: Dx 0 ; Dy 0 ; and Dz 0 .Uniform scaling of the x-, y-, and z-dimensions of the cavity by a factor of s results in the use of fewer calculation cells along the z-direction because the space step is larger along this direction (Dz 0 ¼ s 2 Dz > sDz), requiring a factor of s fewer cells than if the spatial resolution along the z-direction was given by z 0 ¼ sDz.The advantage of fewer cavity cells along the z-direction is not exclusive to explicit FDTD routines that depend on the Courant condition as it could be used to increase the computational efficiency of other prediction methods.To calculate the eigenfrequencies of a scaled rectangular room so that they equal those of the actual room, the speed of sound in air must be scaled using c 0 ¼ sc, such that where f 0 n denotes the eigenfrequencies of the scaled room.
(3) To keep the same characteristic acoustic impedance, the density of the air must be scaled using q 0 o ¼ q o =s.(4) The magnitude of the driving-point mobility for the scaled plate needs to be offset from that of the actual plate.The driving-point mobility, Y dp , of a simply supported isotropic plate is given by where i 2 ¼ À1, S is the surface area of the plate, w 2 p;q ðx; yÞ is the local bending mode shape, and x p;q denotes the angular mode frequency.
Since the mode shapes, eigenfrequencies and loss factors are the same for the actual and scaled plates, the only difference between their transfer mobilities is in the absolute value.Taking the absolute value of the transfer mobility yields Hence the following ratio is expected between the drivingpoint mobility of the scaled and actual plates: This indicates that the magnitude of the scaled driving-point mobility is smaller than that of the actual plate, by a factor of s À4 or by 80 log 10 ðsÞ dB when considering mobility in decibels using 20 log 10 ðjY dp jÞ.Therefore the results should be scaled accordingly, in order to account for this offset in the magnitude of the driving-point mobilities.For the experimental validation in this paper, s ¼ 6 and therefore the shift in level is % 62 dB.One limitation concerns the high-frequency limit for pure bending wave theory.If the thin plate frequency limit for the actual plate is f B % 0:05c L =h, the limit for the scaled plate f 0 B is given by f 0 B ¼ f B =s 2 and the error in the simulation results will increase above this limit. 23

Extension to other plate boundary conditions
For plates with boundary conditions other than simply supported boundaries, it is only necessary to be able to calculate or estimate the eigenfrequencies in order to identify the scaling factor for the z-direction.Eigenfrequencies corresponding to the mode (m, n) of a rectangular thin plate were originally approximated by Warburton 25 for any combination of free (F), simply supported (SS) and clamped (C) boundaries: where the term q mn is given by with constants G x , H x , J x , G y , H y , and J y that are characteristic of each boundary condition involving any permutation of F, C, and SS conditions.These constants only depend on the mode indices, m and n.Hence, when lateral dimensions of the plate are scaled by a factor of s, the eigenfrequencies of the scaled system are given by The term q mn is unaffected by geometric scaling, since G x , G y , H x , H y , J x , and J y only depend on the mode indices m and n and the scaling factors on the terms sL x =sL y cancel out.Hence, to obtain the same eigenfrequencies as the original system, it is necessary to scale the thickness of the plate by a factor of s 2 .It is therefore concluded that the scaling methodology is identical for plates with any combination of free, clamped or simply supported boundaries.

Extension to other topologies
The scaling approach is readily applied to more complex problems involving a number of geometrically parallel thin plates and/or acoustic cavities as illustrated by the examples in Fig. 4.This includes the situation which simulates a sound transmission suite that can be used to determine impact or airborne sound insulation-see Fig. 4(d).When scaling a thin plate, the spatial resolution across its lateral dimensions is scaled by a factor of s, whereas the spatial resolution along the thickness direction must be scaled by a factor of s 2 , so that the eigenfrequencies predicted by Eq. ( 18) remain invariant.This scaling of the spatial resolution can be considered as a coordinate transformation that applies to the entire FDTD model, given by x 0 ¼ sx; y 0 ¼ sy; and z 0 ¼ s 2 z.
For models involving one or more parallel plates, the scaling of their dimensions will be congruent because their lateral dimensions will be scaled by a common factor of s and their thickness direction will be scaled by a common factor of s 2 .Therefore all the scaled plates will preserve the dynamic characteristics of the actual plates.For example, the plates shown in Fig. 4(b) can be simultaneously scaled if the thickness value for each of the two plates is set to h 0 1 ¼ s 2 h 1 and h 0 2 ¼ s 2 h 2 .Additionally, since both plates use the same scaling factor, their lateral dimensions are given by L 0 x1 ¼ sL x1 ; L 0 y1 ¼ sL y1 and L 0 x2 ¼ sL x2 ; L 0 y2 ¼ sL y2 .Future work could consider whether it would be computationally advantageous and feasible to apply the scaling approach to orthogonally arranged plates facing into a cavity, i.e., two plates forming an L-junction.

Numerical efficiency of the scaling approach
As noted in the second bullet point in Sec.II C, the scaling approach requires a factor of s fewer elements than without scaling.In addition to the computational gain from the reduction in the number of elements, there is the additional benefit of being able to use a larger time step.This gain in numerical efficiency can be estimated by assuming that Dx ) Dz and Dy ) Dz such that the corresponding terms for Dx and Dy in the Courant condition [Eq.( 24)] can be omitted.For thin plates this approximation is reasonable because the plate thickness (z-direction) is typically at least one order of magnitude smaller than the lateral dimensions of the plate.Hence the relation between the original time step and the scaled time step is estimated to be Therefore the time step using the scaling approach is larger than that obtained without scaling by a factor of up to s 2 .In conclusion, the scaling approach requires s fewer cells and the relationship between the scaled and original time steps is a maximum of s 2 .Hence the total computational time of the scaled model is estimated to be reduced by a factor up to s Â s 2 ¼ s 3 compared to the computation time needed for the original model.Note that this is the maximum possible reduction in computation time; the actual reduction in computation time will be less than s 3 .The plate is positioned at a height of 0.78 m above the floor using a metal frame-see Fig. 6.The minimum distance between the short edge of the plate and the nearest wall is 0.53 m and between the long edge and the nearest wall is 0.33 m.A simply supported boundary condition is applied around its edges by using a heavy steel frame with pins at 20 mm centres as described by Yin and Hopkins. 24o provide damping similar to the Rayleigh curve, a viscoelastic damping material (Sylomer) is fixed onto the plate surface.By applying different configurations of damping material and measuring the loss factors using the 3 dB downpoints in the magnitude of the driving-point mobility, a diamond-shaped configuration (see Fig. 6) is chosen because the overall damping approximately follows a Rayleigh damping curve below 200 Hz.
Sound pressure measurements inside the room are taken using two grids, one horizontal grid (15 Â 11 positions) and one vertical grid (13 Â 11). Figure 5(a) defines the x-, y-, and z-directions for these grids.The horizontal grid is 0.84 m above the floor and 0.06 m above the plate [Fig.5(b)].The vertical grid is 0.82 m from the back wall and 0.32 m from the edge of the plate [Fig.5(c)].The distances between adjacent grid positions are: 0.2 m in the x-direction for the horizontal grid, 0.18 m in the y-direction for both horizontal and vertical grids, and 0.2 m in the z-direction for the vertical grid, as indicated in Fig. 7.An array of six 1/4 in.microphones (B&K type 4135 with B&K type 2670 pre-amplifiers) is used to measure the sound pressure at the grid points.The maximum uncertainty in each microphone position is estimated to be 61.5 cm in the plane of the horizontal or vertical grid.
A force transducer (B&K type 8200) and accelerometer (B&K type 4393) are connected at the excitation point to monitor the input force.This also allows measurement of the driving-point mobility in order to estimate the modal loss factors.The excitation signal is pseudo-random noise with a cutoff frequency of 200 Hz.
The sound pressure and force signals are used to calculate a complex transfer function of pressure-to-force at all grid points using 0.25 Hz frequency resolution and an upper frequency of 200 Hz.The magnitude of the transfer functions obtained at 0.25 Hz were linearly averaged in order to yield 1 Hz intervals.

B. Numerical implementation
In the FDTD model the material properties for the aluminium plate are: q ¼ 2700 kg/m 3 , c L ¼ 5100 m/s and The properties of air for the acoustic medium were q o ¼ 1.2 kg/m 3 and c ¼ 343 m/s.Using a scaling factor of s ¼ 6 gives q 0 o ¼ 0.20 kg/m 3 and c 0 ¼ 2058 m/s.Previously validated acoustic FDTD models of this room 2 give the specific acoustic impedance for the boundaries as 224.9, which remains the same after scaling.The global Cartesian frame of reference used in the FDTD model is coincident with that shown in Fig. 5.The spatial resolution used for the scaled FDTD model is Dx 0 ¼ 0.39 m, Dy 0 ¼ 0.35 m, and Dz 0 ¼ 0.13 m.The largest wave speed that is accounted for in the vibroacoustic model is the quasi-longitudinal phase velocity of aluminium, 5100 m/s.As discussed in Sec.II C, the largest time step which satisfied the Courant condition and provided stability was found to be 1:93 Â 10 À5 s, which is %15% smaller than the value obtained with Eq. (24).The simulations are carried out over a time interval of 4 s.
The driving function assigned to the source is proportional to the first time derivative of the Gaussian pulse, which has the form where t o is the time offset and r o is the Gaussian width of the pulse and A o is an amplitude constant which was assigned the value of 10 À4 N s 2 =m 2 .This particular waveform is chosen because its spectrum contains no energy at 0 Hz (which would represent static loading).The values chosen for t o and r o determine the frequency content of the source function.In this work, t o ¼ 10 ms and r o ¼ 10 À3 ms. Figure 8 shows the waveform and frequency response of the pulse.It is important that most of the power of the source function lies below the maximum frequency allowed by the domain discretization. 1t is necessary to estimate an upper frequency limit for the FDTD analysis; however, there is more than one factor that determines this limit.The sampling frequency used in the FDTD simulations is 51 724 Hz and according to the Nyquist sampling theorem this results in an upper limit of 25 862 Hz.In contrast, the upper limit due to the scaling approach for the aluminium plate with a scaling factor of 6 is 1418 Hz.In terms of spatial discretization for the air, the use of six computational cells per wavelength results in an upper limit of %870 Hz (based on c 0 ¼ sc ¼ 2040 m/s).For bending waves on the plate, six cells per wavelength gives an upper limit of %300 Hz; hence as this is the lowest value, it provides an estimate for the upper frequency limit of the FDTD simulation.
The FDTD simulation outputs a time history of transient pressure, p(t), and a force driving function, f(t), each with a duration of 4 s and a time step of 1:93 Â 10 À5 s.A fast Fourier transform (FFT) was carried out to give complex PðxÞ and FðxÞ.To obtain the complex transfer function P/F, the complex vectors PðxÞ and FðxÞ were pointwise divided.After obtaining the complex transfer function, its magnitude was taken to give jPðxÞ=FðxÞj from which the level in decibels was calculated using 20 log 10 ðjPðxÞ=FðxÞjÞ.The logarithmic transfer function in 0.25 Hz lines was arithmetically averaged over every four points to give a frequency resolution of 1 Hz.This resolution corresponds to that of the measurements, which was also averaged to 1 Hz from the original data measured at 0.25 Hz.

IV. RESULTS
The FDTD model of the experimental situation uses a scaling factor of s ¼ 6 for which the computation time was reduced by using the scaling approach.After accounting for the total number of calculation cells and the time step used in the scaled model, the computation time was reduced by a factor of %170 compared to the original model.As expected (see Sec. II C 3), this factor does not exceed the ratio of 216 that corresponds to s 3 indicated by Eq. (32).
To assess numerical dispersion in the FDTD model for wave motion in the acoustic medium, a hard, point pressure source is implemented in one corner of the empty room (i.e., without the plate).The sound pressure response in a different corner is then used to identify modal peaks for comparison with the analytical eigenfrequencies that are calculated for an empty room with rigid boundaries.Results for modes below 200 Hz are shown in Table II.These indicate that the errors are less than 2.2%.To assess numerical dispersion for the elastic plate, analytical eigenfrequencies for a simply supported plate are compared with the modal peaks in the FDTD driving-point mobility below 200 Hz.The results are shown in Table III which indicates that the errors are no more than 5.1%.Hence, numerical dispersion can be considered to be negligible for both the air and plate below 200 Hz.
For the plate, a comparison of the frequencies at which the peaks occur in the measured and FDTD driving-point mobilities is shown in Table IV.The agreement confirms that the experimental setup provides a reasonable approximation of simply supported plate boundaries, with the largest difference (7.7%) occurring for mode 2, the f 2;1 mode.
Figure 9 shows the driving-point mobility after accounting for the level offset due to scaling and allows a comparison between FDTD and the measurements.As expected from Table IV, there are slight differences in the frequencies at which the peaks occur but there is reasonable agreement in terms of level with differences ranging from 0.4 to 7 dB.This indicates that the approach used to model the damping of the plate is appropriate, which is also confirmed by comparing loss factors from measurements and FDTD using the 3 dB down points in the driving-point mobility in Table V.The agreement between measured and predicted loss factors is reasonable for the first three modes but has errors of %50% for modes 4 and 5.
To assess the ability of FDTD to predict the spatial variation in sound pressure in the room, a comparison is now made between measured and FDTD transfer functions (magnitudes).The transfer functions for all grid points in the horizontal and vertical grids are shown in Figs. 10 and  At frequencies corresponding to plate modes f 11 and f 12 that occur below the lowest room mode f 010 , the contour plots in Figs. 12 and 13 show close agreement between measurements and FDTD.For the horizontal grid, the sound pressure field corresponds to the vibration field of the plate mode.The results for the vertical grid show that the sound pressure level varies by up to %40 dB over the grid surface.This demonstrates that it is inappropriate to assume a uniform sound field (pressure zone) below the first room mode in a small acoustic cavity which is excited by a plate.
Figure 14 shows the spatial variation above the lowest room mode at a frequency close to the lowest axial mode f 001 (vertical direction) and in between plate modes f 12 and f 21 .In terms of the spatial variation in the horizontal and vertical grids, there is close agreement between measurements and FDTD, with the vertical grid showing the expected variation in sound pressure corresponding to the lowest axial mode in the vertical direction.However, FDTD underestimates the level by %8 dB for both grids.This issue in predicting the correct level has been observed to occur at other frequencies where there is a room mode that is in between plate modes where at least one of the plate modes has a p or q as an even number.This could be due to cancellation in the radiated field that occurs with the unbaffled plate in the FDTD model but does not occur exactly in the experimental setup due to the metal frame that supports the plate.The spatial variation in the horizontal grid is characterised by low sound pressure levels over the surface of the plate because in the vicinity of the plate it prevents the establishment of the mode shape for the lowest axial mode.
Figure 15 shows the response at 82 Hz near the f 12 plate mode which is in between room modes f 001 and f 011 .There is close agreement between measurements and FDTD for the horizontal grid in terms of the spatial variation and levels.However, there is less agreement for the vertical grid, particularly at grid positions that are at a higher elevation than the plate.
Figure 16 shows close agreement between measurements and FDTD for the horizontal and vertical grids for the response at 92 Hz which is close to the lowest tangential room mode f 011 and in between excited plate modes f 21 and f 22 .There is similarly close agreement in Fig. 17  The general finding from the comparison of measured and predicted transfer functions is that FDTD is capable of predicting the spatial variation of sound pressure within experimental error.However, in many noise control situations it is sufficient to predict the spatial-average sound pressure level in the cavity.This also applies to airborne and impact sound insulation involving small rooms at lowfrequencies where the spatial-average corresponds to the average over the entire room volume. 26To indicate the accuracy after spatial-averaging, Fig. 18 shows differences between the measured and FDTD spatial-average magnitude of the transfer functions for all fifteen peaks that occur below 200 Hz.This scatter plot indicates that 60% of the data points are within 63 dB, and that 76% are within 66 dB.
V. CONCLUSIONS Vibroacoustic modelling using FDTD has been developed for an elastic plate undergoing point excitation and radiating into an acoustic cavity.In comparison with room acoustic simulations, it can be computationally expensive to run a large vibroacoustic model with a fine spatial resolution because wavespeeds for structure-borne sound are relatively high.In this paper a scaling approach is proposed and validated to overcome this problem.Modifications to the geometry and physical properties are used to preserve the dynamic characteristics of the model whilst allowing much larger time steps.This reduces the total number of iterations necessary to complete the simulation and significantly reduces computation times.An alternative approach is also proposed and implemented in FDTD to model the boundaries between the air and the solid medium with improved computational efficiency.This results in the model in this paper being %170 times faster than conventional FDTD and the calculations could be parallelized to give further reductions in computation time.The scaling approach can be applied to more complex problems that involve more than one geometrically parallel thin plate and more than one acoustic cavity.
Both modelling approaches proposed in this paper are experimentally validated by the agreement between FDTD and measurements.This confirms the validity of implementing a thin plate undergoing bending wave motion as a threedimensional solid that can support multiple wave types.In the frequency range below the lowest room mode, the close agreement between FDTD and measurements shows the existence of large variations in sound pressure level.This confirms the importance of having a validated vibroacoustic prediction model to predict sound fields inside acoustic cavities in the low-frequency range.

FIG. 2 .
FIG. 2. (Color online) Lattice diagram for a cross-section through the solid medium (shaded grey) indicating the implementation of the simply supported boundary condition using a velocity node set to zero indicated by 1.

FIG. 4 .
FIG. 4. (Color online) Examples of valid configurations for the scaling method: (a) and (b) two isolated, parallel plates, (c) two isolated plates that each face into an acoustic cavity and (d) two acoustic spaces separated by a plate.Light grey surfaces represent the boundaries of the acoustic cavity.

FIG. 11 .
FIG. 11.Transfer functions for all grid points in the vertical grid: (a) measured and (b) FDTD.
11, respectively.Peaks in these transfer functions correspond to global resonances of the plate-cavity system for which there are fifteen peaks below 200 Hz.For the first six of these global resonances, contour plots are shown in Figs.12À17 with the outline of the plate indicated using solid black lines.These plots allow comparisons of measurements and FDTD for the horizontal and vertical grids.
Figure16shows close agreement between measurements and FDTD for the horizontal and vertical grids for the response at 92 Hz which is close to the lowest tangential room mode f 011 and in between excited plate modes f 21 and f 22 .There is similarly close agreement in Fig.17at 104 Hz which is close to the f 22 plate mode and in between room modes f 011 and f 110 .The general finding from the comparison of measured and predicted transfer functions is that FDTD is capable of predicting the spatial variation of sound pressure within experimental error.However, in many noise control

TABLE I .
Ratios of eigenfrequencies obtained using effective and expected boundaries using different numbers of normal stress nodes.

TABLE II .
Empty room: comparison of analytical eigenfrequencies and the frequencies of peaks in the room response from FDTD.

TABLE III .
Plate: comparison of analytical eigenfrequencies and modal peaks in the driving-point mobility from FDTD.

TABLE IV .
Plate: comparison of measured and FDTD eigenfrequencies identified from modal peaks in the driving-point mobility.

TABLE V .
Plate: comparison of measured and FDTD loss factors.
FIG. 10.Transfer functions for all grid points in the horizontal grid: (a) measured and (b) FDTD.