Statistical energy analysis of inhomogeneous systems with slowly varying properties

Journal of Sound and Vibration 333 (2014) 7216–7232

Contents lists available at ScienceDirect

Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi

Statistical energy analysis of inhomogeneous systems with slowly varying properties J. Legault n, J. Woodhouse, R.S. Langley Cambridge University Engineering Department, Trumpington Street, Cambridge CB2 1PZ, UK

a r t i c l e in f o

abstract

Article history: Received 27 April 2014 Received in revised form 7 August 2014 Accepted 19 August 2014 Handling Editor: I. Lopez Arteaga Available online 19 September 2014

This paper is concerned with the analysis of high frequency vibrations of inhomogeneous systems with slowly varying properties subjected to harmonic excitation. A formula for the spatial scaling of the vibrational energy density of such systems is derived and is used to predict the mean and the variance of their vibrational energy response. The derivations are based on the assumption that waves propagate with negligible partial reflections in the systems, which holds when the relative variation of the properties of the systems is small over a wavelength of vibration. To validate the theoretical findings, numerical simulations of two plate systems are performed: a plate with a linearly varying thickness and a plate with a thickness bump. & 2014 Elsevier Ltd. All rights reserved.

1. Introduction The analysis of high frequency vibrations of built-up engineering structures such as cars or aircraft is a challenging task. Standard analysis techniques such as the finite element method require a large number of degrees-of-freedom (sometimes millions) to describe the vibration response of the structure accurately. Moreover, the vibration response becomes sensitive to variability in the structure's properties (geometry, material characteristics, boundary conditions, etc.), and thus predicting the vibration response of a structure with fixed properties is of little use. Instead, the analyst must predict the statistics of the vibration response of an ensemble of structures whose properties will vary slightly from the properties of the nominal structure. To tackle this challenge, energy-based methods are often used. The most prevalent of these methods is Statistical Energy Analysis (SEA) [1]. In SEA, the structure is divided into subsystems, and the ensemble statistics of the total timeaveraged vibrational energy in each subsystem are sought. An important task when building an SEA model is to define the subsystems. In order to yield accurate predictions, the subsystems must be weakly coupled. Unless the damping is very small, weak coupling generally occurs when the subsystems are weakly connected, i.e. when energy cannot flow freely between the subsystems because there is a significant impedance mismatch at their interface [2,3]. Defining the subsystems is often straightforward (e.g. two plates with different thicknesses coupled along an edge), but some systems are more problematic. A particularly awkward situation is when the properties of the system vary slowly in space, in the sense that their relative variation is small over a wavelength of vibration. Beams with slowly varying cross-section or plates with slowly varying thickness or curvature are good examples. The fact that the properties vary slowly is problematic because the waves will propagate with negligible partial reflections and energy will flow freely. Therefore, according to conventional wisdom, inhomogeneous systems with slowly varying

n

Corresponding author. Present address: Wave6 LLC, 11965 El Camino Real, San Diego, CA 92122, USA. E-mail address: [email protected] (J. Legault).

http://dx.doi.org/10.1016/j.jsv.2014.08.026 0022-460X/& 2014 Elsevier Ltd. All rights reserved.

J. Legault et al. / Journal of Sound and Vibration 333 (2014) 7216–7232

7217

properties cannot be split into several SEA subsystems since the coupling between the subsystems will be strong. Instead, these systems should be treated as single SEA subsystems. There is a vast and rich literature about wave propagation in inhomogeneous media, and the vibration response of inhomogeneous systems with slowly varying properties has been studied for a wide range of systems such as beams and plates [4], shells [5], acoustic wave-guides modeling the energy flow in the cochlea [6], or power law profiled plates exploiting the acoustic “black hole” effect to increase damping [7]. Most often, the WKB approximation (after Wentzel, Kramers and Brillouin) was used for the analysis. However, despite all this research, the statistical energy analysis of subsystems with slowly varying properties remains an open problem. More precisely, closed-form expressions that predict how the ensemble average of the vibrational energy density and the mode shapes of such systems scale with their spatially varying properties have yet to be established (in homogeneous systems, these quantities are assumed to be uniform, although response concentration effects can occur near the boundaries of the subsystem or near the driving location when the subsystem is excited by a localized load [1]). The objective of this paper is to address this shortcoming. It should be noted that several methods besides SEA exist for the prediction of the vibrational energy density of complex structures [8], and some of these methods such as Statistical modal Energy distribution Analysis [9] or Dynamical Energy Analysis [10] could be used to predict the scaling of the vibrational energy density in inhomogeneous subsystems with slowly varying properties. However, these methods would yield numerical predictions, not closed-form expressions. By seeking to derive such expressions, the present paper is thus within the spirit of traditional SEA. The paper is structured in two main parts: the theoretical development (Section 2) and the numerical validation (Section 3). In Section 2, a proportionality relation for the scaling of the vibrational energy density with the properties of the system is obtained by using two different arguments: a subsystem-based argument and a wave-based argument. It is shown that both arguments lead to the same relation and that this relation is consistent with existing results in the literature obtained from the WKB approximation. Formulas for the statistics of the mode shapes, of the total vibrational energy response and of the point-to-point transfer function of the system are then derived. In Section 3, these theoretical findings are validated by studying the energy response of two plate systems: a plate with a linearly varying thickness and a plate with a thickness bump.

2. Theoretical development 2.1. Preliminary considerations The following analysis is concerned with the vibration response of a spatially inhomogeneous, time-invariant, linear dynamical system. The properties of the system are assumed to vary slowly in space, in the sense that their relative variation is small over a wavelength of vibration. The system has a scalar response variable wðx; tÞ representing the displacement of the system at spatial location x and time t. The dimension of x corresponds to the physical dimension of the system, so that, for example, a plate would be described by the coordinates x ¼ ðx1 ; x2 Þ. The system is taken to have light structural proportional damping, with loss factor η, and a distributed harmonic load of frequency ω and complex amplitude PðxÞ is applied. The steady-state response of the system can be expressed as a modal sum in the form [11]

 P r ϕr ðxÞ iωt wðx; t Þ ¼ RefW ðxÞeiωt g ¼ Re ∑U r ϕr ðxÞeiωt ¼ Re ∑ 2 e ; (1) 2 r r ωr ð1 þiηÞ  ω where WðxÞ is the complex displacement amplitude and Ur, ϕr ðxÞ, ωr and Pr are respectively the complex mode amplitude, mode shape, natural frequency and modal load associated with the rth mode of vibration. The modal load Pr is given by Z (2) P r ¼ PðxÞϕr ðxÞ dx; S

where the integration region S corresponds to the spatial domain of the system (length, area or volume). The mode shapes are taken to be scaled to unit generalized mass so that
Z 0 if r a r 0 μðxÞϕr ðxÞϕr0 ðxÞ dx ¼ ; (3) 1 if r ¼ r 0 S where μðxÞ is the appropriate mass density of the system (for example, the mass per unit area for a plate or the mass per unit length for a beam). The total time-averaged vibrational energy of the system, which is equal to twice the total timeaveraged kinetic energy for a resonant system, can be derived from Eqs. (1)–(3) in the form Z E¼

a r ω2 eðxÞ dx ¼ 2∑ ; 2 2 r ðωr  ω2 Þ þ η2 ω4 S r

(4)

1 2 ω μðxÞjW ðxÞj2 2

(5)

where eðxÞ ¼

7218

J. Legault et al. / Journal of Sound and Vibration 333 (2014) 7216–7232

is the time-averaged vibrational energy density and ar ¼

1 jP r j2 4

(6)

the modal coefficient. From here onwards, the total time-averaged vibrational energy and the time-averaged vibrational energy density will be simply called the total energy and the energy density. This study is concerned with the statistics of the response of an ensemble of random systems whose properties vary slightly from that of the nominal inhomogeneous system, so that the natural frequencies and the mode shapes which appear in Eqs. (1)–(6) are taken to be random variables. Throughout the paper, E½ is used to represent averaging over this ensemble. To simplify the terminology, the ensemble average of a given quantity X will be simply called the “mean” of X. For a spatially homogeneous system such as a uniform plate or an acoustic volume, it is well known that if the system is sufficiently random then the statistical distribution of its mode shapes and natural frequencies conform to those of the Gaussian Orthogonal Ensemble (GOE) [12–14], regardless of the details of the underlying randomness. In this paper, it is shown that this remains valid for an inhomogeneous system with slowly varying properties. The only difference is that the 2 mean squared value of the mode shapes, i.e. the mean of the quantity ϕr ðxÞ, which will here be called the modal norm, scales with the local properties of the system instead of being spatially uniform. The formula for this scaling is obtained in Section 2.3. It is based on the expression derived in the following subsection for the scaling of the mean energy density. 2.2. Scaling of the mean energy density 2.2.1. Subsystem-based argument Suppose that the inhomogeneous system described in the previous subsection is split into N SEA subsystems with arbitrary boundaries. As mentioned in the Introduction this would normally be avoided since the subsystems will be strongly coupled. However there is a slight twist to this argument, as follows. Let E1 ; E2 ; …; EN be the total energies of the N subsystems and Π in;1 ; Π in;2 ; …; Π in;N be the powers injected into these subsystems by the loading PðxÞ. In an SEA model, these quantities are related by the following set of N equations [1]:   N     E½Ej  E½Ek  ωηE Ej þ ωnj ∑ ηjk  ¼ E Π in;j ; j ¼ 1; 2; …; N; (7) nk nj kaj where nj is the modal density of subsystem j and ηjk is the coupling loss factor between subsystems j and k. The coupling loss factors satisfy the consistency relationship ηjk nj ¼ ηkj nk , and the modal energy Ej =nj can be written in the form Ej Sj ¼ 〈e 〉; nj nj j

(8)

where 〈ej 〉 and Sj are the spatially averaged energy density and the spatial domain (length/area/volume) of subsystem j, respectively. For the system under consideration, the damping η is small and the subsystems are very strongly connected (waves propagate with negligible partial reflections). The subsystems are thus very strongly coupled. In the strong coupling limit, i.e. when ηjk ⪢η 8 j; k, Eq. (7) predicts equipartition of the mean modal energies of the subsystems [15]: E½〈e1 〉 E½〈e2 〉 E½〈eN 〉 ¼ ¼… : n1 =S1 n2 =S2 nN =SN

(9)

Although SEA does not yield accurate predictions when the coupling is moderately strong, Eq. (9) will be accurate in the limit. The inaccurate predictions are obtained only in the transition between weak coupling and equipartition. Now, consider the situation when the number of subsystems is increased significantly such that their size becomes very small compared to the wavelength at frequency ω. In that situation the SEA model is no longer applicable since the subsystems bear no modes of vibration in the vicinity of ω. However if ω is increased such that the wavelength remains an order of magnitude smaller than the size of the subsystems, the SEA model remains valid. Hence, in the high subsystem and high frequency limit (N-1 and ω-1), Eq. (9) leads to the following asymptotic result: ~ E½eðxÞ ¼ AnðxÞ;

(10)

~ where nðxÞ is the local modal density and A is a positive constant that does not depend on x. The mean of the local modal ~ energy eðxÞ=nðxÞ is thus expected to be spatially uniform in the high frequency limit. ~ The local modal density nðxÞ can be obtained by calculating the modal density of a homogeneous system that has the properties of the inhomogeneous system at location x, normalized per unit length/area/volume. For example, the local modal density of an inhomogeneous (but locally isotropic) thin plate in bending with variable thickness hðxÞ is given by [11] !1=2   1 μðxÞ 1=2 1 12ð1  ν2 Þρ ~ n ðx Þ ¼ ¼ ; (11) 2 4π DðxÞ 4π Eh ðxÞ where DðxÞ, ρ, E and respectively.

ν are bending stiffness, mass density per unit volume, Young's modulus and Poisson ratio of the plate

J. Legault et al. / Journal of Sound and Vibration 333 (2014) 7216–7232

7219

To obtain the global modal density n of the system, the local modal density can be integrated over its spatial domain: Z ~ dx: (12) n ¼ nðxÞ S

Eq. (12) follows directly from the semi-classical formula for the mean density of energy states of a quantum mechanical system [16], which effectively integrates the local modal density over the system domain. Combining Eqs. (4), (10) and (12) leads to the following relation between the energy density and the total energy: E½eðxÞ ¼

~ nðxÞ E½E: n

(13)

Eq. (13) shows that for a system with slowly varying properties, the energy density simply scales with the ratio of the ~ local modal density to the global modal density. If the properties of the system are spatially uniform, nðxÞ ¼ n=S and Eq. (13) yields E½eðxÞ ¼ E½E=S, which is the standard SEA result for a homogeneous system.

2.2.2. Wave-based argument The following analysis presents an alternative argument to the subsystem-based argument described in the previous subsection. The argument is based on energy flow conservation in a diffuse wave field. The derivation is presented for an arbitrary two-dimensional inhomogeneous system, but it may be extended to three-dimensional systems or simplified to one-dimensional systems without loss of generality. The dispersion equation of a two-dimensional inhomogeneous system can be established by considering the propagation of a plane wave with the space–time dependency eiωt  ik1 x1  ik2 x2 , where k1 and k2 are the wavenumbers in the x1 and x2 directions. By substituting this form of dependency in the differential equation of motion of the system, and by neglecting the terms associated with the spatial variation of the properties of the system, a relation of the form

ω ¼ f ðk; xÞ ¼ f ðk1 ; k2 ; x1 ; x2 Þ

(14)

is obtained. The terms associated with the spatial variation of the properties of the system are neglected because their relative variation is assumed to be small over a wavelength of vibration. Eq. (14) is equivalent to the Eikonal equation in the WKB approximation. In the wavenumber plane, the ðk1 ; k2 Þ contour line obtained by fixing ω and x in the dispersion equation is called the local dispersion curve. This curve can be parameterized by switching to polar coordinates, i.e. by making the change of variable k1 ¼ k cos θ and k2 ¼ k sin θ, where θ is the wave heading angle and k is the (heading dependent) scalar wavenumber. As shown in Fig. 1a, the local dispersion curve can also be parameterized by using the local line coordinate s. Hence, the scalar wavenumber depends on both x and s in the local dispersion equation. It should be noted that for conciseness, only the first quadrant of the wavenumber plane is shown in Fig. 1a. To properly display the local dispersion curve of an arbitrary anisotropic system, at least two quadrants are required since the dispersion curve has only one axis of symmetry in the wavenumber plane [17]. The phase velocity c and the group velocity vector cg associated with the wave component at coordinate s are given by cðs; xÞ ¼ ω=kðs; xÞ and   ∂ω ∂ω cg ðs; xÞ ¼ ∇ω ¼ ; ; (15) ∂k1 ∂k2 respectively. By definition the group velocity vector is perpendicular to the dispersion curve and the amplitude of its

Fig. 1. Dispersion curve in the wavenumber plane. (a) At location x. (b) At locations xa and xb .

7220

J. Legault et al. / Journal of Sound and Vibration 333 (2014) 7216–7232

projection on the heading direction is given by cgθ ðs; xÞ ¼ cg ðs; xÞ cos ðθg  θÞ ¼ ∂ω=∂k; where cg ðs; xÞ ¼ ‖cg ðs; xÞ‖ is the amplitude of the group velocity vector, and Following Langley [17], the local modal density of the system is given by Z π 1 kðθ; xÞ n~ ðxÞ ¼ 2 dθ : 2π 0 cgθ ðθ; xÞ

(16)

θg is the angle between cg and the k1-axis. (17)

When the loading PðxÞ is applied to the system, a wave field that satisfies the local dispersion equation at every point will be generated (waves will adopt the local wavenumber as they propagate), and the wave component ðk1 ; k2 Þ associated with ^ xÞ. In SEA, it is assumed that the wave field is diffuse, line coordinate s will have an energy density per unit s equal to eðs; which means that different wave components contribute independently to the mean energy [18]: Z ^ xÞ ds; E½eðxÞ ¼ E½eðs; (18) Γ

where Γ corresponds to the dispersion curve contour in the wavenumber plane. For spatially homogeneous isotropic ^ xÞ ¼ e^ 0 , systems, a diffuse field also implies that waves have the same energy in all directions and at all locations, i.e. E½eðs; ^ xÞ is where e^ 0 is a constant that does not depend on x or s. In the case of spatially homogeneous anisotropic systems, E½eðs; independent of x but it is not independent of s as shown by Langley [19]. The objective of the following analysis is to extend the mathematical definition of a diffuse field to spatially inhomogeneous anisotropic systems and to study the implications. Consider two arbitrary locations xa and xb in an inhomogeneous anisotropic system. In the wavenumber plane the dispersion curves at these two locations can be described by using the local line coordinates sa and sb as shown in Fig. 1b. The corresponding group velocity vectors are denoted cg ðsa ; xa ) and cg ðsb ; xb Þ. These vectors make angles θg;a and θg;b with the k1-axis. The small dispersion curve segments associated with the infinitesimal segment dk1 centred on k1 are noted dsa and dsb . Since the partial reflections of the waves in the system are assumed to be negligible, the wave field generated by the external loading can be assumed to be diffuse, and hence the mean flow of wave energy has to be conserved for each wave component. Therefore the mean flow of wave energy of the wave component for fixed k1, projected in the direction k2, should be the same at locations xa and xb . Mathematically this condition can be written in the form [18] ^ b ; xb Þcg ðsb ; xb Þ sin ðθg;b Þ dsb : ^ a ; xa Þcg ðsa ; xa Þ sin ðθg;a Þ dsa ¼ E½eðs E½eðs

(19)

Noting that sin θa dsa ¼ sin θb dsb ¼ dk1 , Eq. (19) yields ^ b ; xb Þcg ðsb ; xb Þ: ^ a ; xa Þcg ðsa ; xa Þ ¼ E½eðs E½eðs

(20)

For Eq. (20) to be true for any arbitrary locations xa and xb , the following condition must thus be satisfied at every location x: ^ xÞcg ðs; xÞ ¼ B; E½eðs;

(21)

where B is a positive constant that does not depend on s or x. Combining Eqs. (18) and (21), switching to polar coordinates and using Eqs. (16) and (17) yield Z 2π Z 1 kðθ; xÞ ds ¼ B E½eðxÞ ¼ B dθ ¼ An~ ðxÞ; (22) cg ðθ; xÞ cos ðθg  θÞ Γ cg ðs; xÞ 0 where A is a positive constant that does not depend on x. Similar to the SEA subsystem-based argument (see Eq. (10)), equipartition of local modal energy is obtained. It should be noted that the proportionality relation of Eq. (21) is invariant to the choice of the axes k1 and k2, in the sense that the same result would be obtained by performing the analysis for fixed k2 instead of fixed k1, or by rotating the axes k1 and k2 by an arbitrary angle in the wavenumber plane. 2.2.3. Link with the WKB approximation To the best knowledge of the authors, a scaling relation that has the form of Eq. (22) has not been presented in the literature. However, this relation can be related to the transport equation obtained when applying the WKB approximation, and it is well known that the transport equation governs the high-frequency energy density of the wave field, and that this equation is consistent with the principle of conservation of energy flow. More precisely, the transport equation predicts that the energy flow will be conserved along a narrow tube of rays [5,20]. For two- and three-dimensional systems, this conservation law is different from that of Eq. (22); it expresses conservation of energy flow for a single system in a direction that varies as the rays are bent due to the spatially varying properties of the system. In contrast, Eq. (22) expresses conservation of mean energy flow for a fixed wave component in a fixed direction over an ensemble of systems. For onedimensional systems however, the two conservation laws should be consistent since there is only one direction for the traveling waves. It can be shown that this is indeed the case. For instance, Pierce [4] applied the WKB approximation to an Euler–Bernoulli beam with slowly varying properties and showed that for the flow of wave energy to be conserved, the amplitude jWðxÞj of the waves must scale as jWðxÞj ¼ Q ½μðxÞ  3=8 ½DðxÞ  1=8 ;

(23)

J. Legault et al. / Journal of Sound and Vibration 333 (2014) 7216–7232

7221

where Q is a constant that does not depend on x, and where μðxÞ and D(x) are the mass per unit length and the bending stiffness of the beam, respectively. Since the local modal density of an Euler–Bernoulli beam is given by [11]   1 μðxÞ 1=4 n~ ðxÞ ¼ ; (24) 2πω1=2 DðxÞ combining Eqs. (5), (23) and (24) yields ~ eðxÞ ¼ πω5=2 Q 2 nðxÞ;

(25)

which is consistent with Eq. (22). 2.3. Mean modal norm In Section 2.2.2, a wave point of view was adopted to determine how the energy density scales with the properties of the system when the response wave field is assumed to be diffuse. If a modal point of view is instead adopted, a diffuse field implies that for resonant modes, (i) the modal amplitudes are uncorrelated zero-mean random variables of equal energy, (ii) the mode shapes are zero-mean random variables of equal variance and (iii) the modal amplitude and the mode shape of a single mode are uncorrelated [21]. Under such conditions, Eqs. (1) and (5) yield E½eðxÞ ¼ I μðxÞE½ϕr ðxÞ; 2

(26)

where I is a positive constant that does not depend on x. Combining Eqs. (3), (4) and (26) gives I ¼ E½E, and therefore it follows from Eqs. (13) and (26) that E½ϕr ðxÞ ¼ 2

~ nðxÞ : nμðxÞ

(27)

The mean of the modal norm ϕr ðxÞ is found to be proportional to the local modal density and inversely proportional to 2 the local mass density. For a homogeneous system, the standard result E½ϕr ðxÞ ¼ 1=M, where M corresponds to the total mass of the system, is recovered. Assuming that GOE statistics also apply to the mode shapes of inhomogeneous systems with slowly varying properties (this issue will be investigated numerically in Section 3.3), the mode shapes are thus expected to be zero-mean uncorrelated Gaussian random variables [12,14] with mean squared value given by Eq. (27) and mode shape factor KðxÞ given by 4 2 KðxÞ ¼ E½ϕr ðxÞ=E2 ½ϕr ðxÞ ¼ 3. 2

2.4. Energy response statistics Having derived a closed-form expression for the mean modal norm, the statistics of the energy response of the system can now be predicted. Under GOE statistics, point-process theory leads to the following expressions for the ensemble mean and relative variance of the total energy [14]:

πn ; ηω

(28)

α þ qðmÞ ; πm

(29)

E½E ¼ E½ar  relvar½E ¼

where m ¼ ηωn is the (global) modal overlap factor, α is the loading spatial factor and q(m) is a function of the modal overlap factor. This function is given by Eq. (24) of Ref. [22]. It tends towards zero at low values of m and towards  1 at high values of m. If a point force of complex amplitude F0 is applied at location x0 , ¼ E½a2r =E2 ½ar 

E½ar  ¼

jF 0 j2 2 E½ϕr ðx0 Þ; 4

(30)

where E½ϕ is given by Eq. (27). Moreover, α ¼ Kðx0 Þ ¼ 3 (see Section 2.3). Therefore, the mean of the total energy will depend on the location of the excitation but the relative variance will not. This conclusion is important since it could have ~ ~ been thought that the relative variance depends on the local modal overlap factor defined by mðxÞ ¼ ηωnðxÞ. A quantity that is closely related to the total energy under point forcing and that is important for power input considerations is the real part of the input mobility YðxÞ, which is defined as the ratio of the velocity to the force at the drive point. Indeed, it is known that for point forcing [1], 2 r ðx 0 Þ

E½E ¼

jF 0 j2 E½RefY ðxÞg: 2ηω

(31)

Combining Eqs. (27), (28), (30) and (31) leads to E½RefY ðxÞg ¼

~ π nðxÞ : 2μðxÞ

(32)

7222

J. Legault et al. / Journal of Sound and Vibration 333 (2014) 7216–7232

For a homogeneous system, the standard SEA result E½RefYðxÞg ¼ π n=2M [23] is recovered. Hence, Eq. (32) shows that the mean of the real part of the input mobility of an inhomogeneous system with slowly varying properties corresponds to the local input mobility, i.e. to the mobility of a homogeneous system that has the local properties of the inhomogeneous system. Although intuitive, this conclusion is important since input mobility measurements are frequently used in experimental SEA. The other type of loading that is considered in this study is rain-on-the-roof forcing. This type of loading is known to generate diffuse fields [24] and is thus commonly assumed in SEA. It is a spatially delta-correlated random loading such that 〈PðxÞP n ðx0 Þ〉P ¼

S^ pp μðxÞδðx x0 Þ; M

(33)

where S^ pp is the squared magnitude of the forcing spectrum and where 〈  〉P represents averaging over all loading realizations. Combining Eqs. (1)–(6) and (33) leads to the following expressions for the modal coefficient and the squared magnitude of the displacement amplitude under such loading: ar ¼

jW ðxÞj2 ¼

S^ pp ; 4M

S^ pp ϕ2r ðxÞ ∑ : M r ðω2r  ω2 Þ2 þ η2 ω4r

(34)

(35)

In the above expressions, the ensemble average over all loading realizations is considered implicit. To simplify the terminology, the squared magnitude of the displacement amplitude will be simply called the “squared displacement” from here onwards. If the total energy is averaged over all loading realizations prior to the calculations of the ensemble mean and relative variance, it follows from Eq. (34) that E½ar  ¼ S^ pp =4M and α ¼1 in Eqs. (28) and (29). Moreover, by comparing Eq. (35) with Eq. (4), it can be seen that the squared displacement under rain-on-the-roof forcing is the reciprocal case of the total energy under point forcing, in the sense that there is only a deterministic factor between the two quantities. Therefore, it follows from Eqs. (28) and (35) that the mean squared displacement under rain-on-the-roof forcing is given by   S^ pp π n 2 E jW ðxÞj2 ¼ E½ϕr ðx0 Þ; M 2ηω3

(36)

while the relative variance can be obtained by replacing the energy with the squared displacement in Eq. (29) and by using α ¼ KðxÞ ¼ 3. The statistics of the squared displacement under point forcing are considered in the next subsection. 2.5. Transfer function statistics By definition, if a point force of complex amplitude F0 is applied at location x0 , the mean squared displacement at location x will be given by E½jWðxÞj2  ¼ jF 0 j2 E½jHðx; x0 Þj2 , where

ϕr ðxÞϕr ðx0 Þ ω2r ð1 þ iηÞ  ω2

Hðx; x0 Þ ¼ ∑ r

(37)

is the point-to-point transfer function of the system, i.e. Green's function of the system in the frequency domain. Adopting a wave point of view, the transfer function can also be written as the sum of a direct field component and a reverberant field component: Hðx; x0 Þ ¼ H dir ðx; x0 Þ þH rev ðx; x0 Þ:

(38)

The reverberant field component corresponds to the contribution from the multiple wave reflections from the boundaries while the direct field component corresponds to the direct field (or free field) Green's function, i.e. to Green's function if the system was infinite or had perfectly absorbing boundaries. Shorter and Langley [25] have shown that when there is sufficient uncertainty in the system boundary conditions, the direct and reverberant field components are uncorrelated, and the reverberant component averages to zero, i.e. E½Hðx; x0 Þ ¼ H dir ðx; x0 Þ:

(39)

The mean squared transfer function magnitude can thus be written in the form E½jHðx; x0 Þj2  ¼ jH dir ðx; x0 Þj2 þ E½jH rev ðx; x0 Þj2 :

(40)

Assuming GOE statistics, Langley and Cotoni [22] derived closed form expressions for the direct field component and for the mean squared magnitude of the reverberant field component. Recently, Choi et al. [26] extended Langley and Cotoni's results by taking the real part of the direct field component into account. Choi et al. also considered the influence of coherent reflections at the boundaries and of randomly distributed scatterers. These effects are not taken into account here, but when they are, Choi et al. showed that E½Hðx; x0 Þ a H dir ðx; x0 Þ. In both Refs. [22] and [26], the expressions were derived for a homogeneous isotropic system. The objective of the following analysis is to extend the results of Refs. [22] and [26] to

J. Legault et al. / Journal of Sound and Vibration 333 (2014) 7216–7232

7223

inhomogeneous systems with slowly varying properties. The analysis is presented for a two-dimensional thin plate in bending but it can be generalized to other types of 2D systems or to 3D systems. As shown in Refs. [22] and [26], the prediction of H dir ðx; x0 Þ and E½jH rev ðx; x0 Þj2  involves the correlation function of the mode shapes, defined as C ðx; x0 Þ ¼

E½ϕr ðxÞϕr ðx0 Þ ðE½ϕr ðxÞE½ϕr ðx0 ÞÞ1=2 2

2

:

(41)

It follows from wave-duality arguments (see e.g. [16]) that the appropriate value of Cðx; x0 Þ is the diffuse field correlation function, which for homogeneous 2D systems has the form [23] Cðx; x0 Þ ¼ CðkRÞ ¼ J 0 ðkRÞ;

(42)

where J is the Bessel function of the first kind, R ¼ ‖x  x0 ‖ is the distance between x and x0 , and k ¼ ðω2 μ=DÞ is the bending wavenumber in the plate, evaluated at ω ¼ ωr . For inhomogeneous systems with slowly varying properties, a closed-form expression like Eq. (42) has not been found in the literature. However progress can be made by noting that the diffuse field correlation function is known to be proportional to the imaginary part of the direct field Green's function [24], i.e. 1=4

Cðx; x0 Þ ¼ A0 ImfH dir ðx; x0 Þg;

(43)

0

where A is a real constant that does not depend on x or x0 . For a homogeneous isotropic thin plate in bending, H dir ðx; x0 Þ can be written in the form [26] H dir ðx; x0 Þ ¼ 

iπ n iπ n 2 2 ðE½ϕr ðxÞE½ϕr ðx0 ÞÞ1=2 Ψ ðkRÞ ¼  Ψ ðkRÞ; 2ω 2ωM

(44)

where n ¼ ðS=4π Þðμ=DÞ1=2 is the global modal density of the homogeneous plate and where Ψ ðkRÞ represents an axisymmetric cylindrical generated at the drive point. This wave function is given by

Ψ ðkRÞ ¼ H0ð2Þ ðkRÞ  Hð2Þ 0 ð  ikRÞ;

(45)

ð2Þ

where H ðÞ is the Hankel function of the second kind. For an inhomogeneous plate with slowly varying properties, H dir ðx; x0 Þ will have a more complicated form since the axisymmetric waves generated at the drive point will be distorted (waves will be bent as they travel). Therefore H dir ðx; x0 Þ will not be axisymmetric. For instance, when the waves approach a region in the system that has a higher stiffness to mass density ratio, the wavefronts will progress faster within this region. Yet, as mentioned in Section 2.2.2, it is expected that waves in such systems will adopt the local wavenumber as they propagate, and that partial reflections will be negligible. Therefore, it is conjectured here that Eq. (44) remains valid, at least approximately, for inhomogeneous plates with slowly varying properties, with the difference that the phase change argument kR is naively replaced with the corrected phase change Z fkRgcorr ¼ kðRÞ dR; (46) R0

where R0 is the straight path between locations x and x0 , and k(R) corresponds to the spatially varying wavenumber along 2 that path. Moreover, the mean modal norm E½ϕr ðxÞ, which is equal to 1/M for a homogeneous plate, is replaced with the prediction of Eq. (27). Making these substitutions in Eq. (44) yields   ~ nðx ~ 0 Þ 1=2 iπ nðxÞ H dir ðx; x0 Þ ¼  Ψ ðfkRgcorr Þ; (47) 2ω μðxÞμðx0 Þ ~ where the local modal density nðxÞ is given by Eq. (11). It thus follows from Eqs. (41), (43) and (47) that     1 nðxÞ ~ nðx ~ 0 Þ 1=2 J 0 ðfkRgcorr Þ; E ϕr ðxÞϕr ðx0 Þ ¼ n μðxÞμðx0 Þ

(48)

where the global modal density n is given by Eq. (12). The phase conjecture described above is based on the assumption that the contour levels of H dir ðx; x0 Þ, i.e. the wavefronts generated at the driving location, can be traced with the corrected phase prediction of Eq. (46); the wavefronts simply correspond to the spatial contour levels for constant phase. The conjecture also assumes that the decay of the amplitude of the real and imaginary parts along the contour levels will be correctly captured by the wave function Ψ. If the wavefronts are not too distorted and the distance between locations x and x0 is not too large compared to the wavelength, the conjecture is expected to work better. This issue will be investigated numerically in Section 3.6. Having obtained an expression for the correlation of the mode shapes, the mean squared magnitude of the reverberant field component can now be derived. Following Eq. (48) and Section 5.3 of Ref. [26], it can be shown that   E jH rev ðx; x0 Þj2 ¼

~ nðx ~ 0Þ  π nðxÞ 1 þjΨ ðfkRgcorr Þj2 ½2 þ qðmÞ : 2ηω3 nμðxÞμðx0 Þ

(49)

7224

J. Legault et al. / Journal of Sound and Vibration 333 (2014) 7216–7232

In the far field, i.e. when the distance between points x and x0 becomes large, jΨ j-0, and the mean squared magnitude of the transfer function simplifies to E½jHðx; x0 Þj2 far ¼

~ nðx ~ 0Þ π nðxÞ : 2ηω3 nμðxÞμðx0 Þ

(50)

~ ~ 0 Þ=μðx0 Þ ¼ n=M, and Eq. (15) of Ref. [22] is recovered from Eq. (50). For a homogeneous plate, nðxÞ= μðxÞ ¼ nðx 3. Numerical validation 3.1. Simulation parameters To validate the theoretical findings, numerical simulations of the bending vibrations of an inhomogeneous (but locally isotropic) rectangular thin plate with simply supported edges have been performed. The governing differential equation of motion of the plate is given by [27] ∂D ∂ð∇2 wÞ ∂D ∂ð∇2 wÞ þ2 ∂x1 ∂x1 ∂x2 ∂x2 ! ! ∂2 w ∂2 D ∂2 D ∂2 w ∂2 D ∂2 D þ 2 þν 2 þ 2 þν 2 ∂x1 ∂x21 ∂x2 ∂x2 ∂x22 ∂x1 D∇4 w þ2

∂2 w ∂2 D ∂2 w þ μ 2 ¼ P; þ2ð1  νÞ ∂x1 ∂x2 ∂x1 ∂x2 ∂t

(51)

where w is the out-of-plane displacement, ν is the Poisson ratio, P is the applied loading, μ is the mass density and D is the bending stiffness. In Eq. (51) it is implicit that these quantities depend on x1 and x2 except the Poisson ratio that is assumed to be uniform. For the simulations, two cases are considered: a plate with a linearly varying thickness (tapered plate) and a plate with a Gaussian thickness bump (bumped plate). These cases are representative of inhomogeneous systems that can be found practice. Beams and plates with slowly varying cross-section or thickness can be found in a wide range of built-up structures while plates with local thickness bumps are often used to provide strength for structural attachments. The tapered and bumped plates are taken to be variations from a spatially uniform rectangular isotropic plate with side lengths L1 and L2, thickness h0, Young's modulus E0, Poisson ratio ν0 and mass density per unit volume ρ0. Therefore, the properties of the tapered plate are the same as those of the uniform plate except the thickness h that varies linearly along the x1-axis as follows:  

 κ  1 2x1 1 : (52) hðxÞ ¼ h0 1 þ κ þ1 L1 In Eq. (52), κ is the taper ratio of plate, i.e. the ratio of the maximal to the minimal thickness. Eq. (52) is such that the total mass M of the tapered plate remains equal to the total mass M0 of the uniform plate regardless of the taper ratio κ. For the bumped plate, the thickness h varies as follows: h γ  ððx1  L1 =2Þ2 þ ðx2  L2 =2Þ2 Þ=2σ2 L1 L2 i hðxÞ ¼ h0 1  γ þ ; (53) e 2πσ 2 where

σ2 ¼

γ

; 2π ðβ 1  γ Þ

(54)

3

Normalized thickness

Normalized thickness

and where γ and β are the bump mass and bump peak ratios respectively. The bump mass ratio is the ratio of the mass contained in the bump to the total mass of the plate while the bump peak ratio is the ratio of the peak thickness of the bump

2 1 0 1

1 0.5

x (m 2 )

0

0

0.5 ) x (m 1

3 2 1 0 1

1 0.5 x (m) 2

0.5 0

0

Fig. 2. Normalized thickness: (a) Tapered plate; (b) bumped plate.

x 1 (m)

J. Legault et al. / Journal of Sound and Vibration 333 (2014) 7216–7232

7225

to the average thickness of the plate. For γ o0:3 and β 42, the bumped plate total mass M is virtually equal to M0 (less than 0.1 percent difference). For the simulations, the following numerical values are used: L1 ¼ 1 m, L2 ¼ 1:2 m, h0 ¼ 0:002 m, E0 ¼ 70 GPa, ν0 ¼ 0:33, ρ0 ¼ 2700 kg m  3 , κ ¼4, γ ¼0.2 and β ¼3. The normalized thickness hðxÞ=h0 of the plates obtained with these values is plotted in Fig. 2. The maximal frequency of interest is 1000 Hz and the damping η is set to 0.01. The modal density n0 of the uniform plate is equal to 0.0307 (rad/s)  1. In the simulations, the plates are randomized with 40 randomly placed point masses totalling 15 percent of the uniform plate mass and with 5 randomly placed rotational springs of stiffness 1000 N m/ rad at each edge. An ensemble containing 5000 realizations has been generated to construct the response statistics, with the location of the point masses and the edge springs being selected from uniform distributions. Throughout the following subsections, the results obtained from this ensemble of random plates are referred to as the “empirical” results. The stiffness and mass matrices of the plates are obtained by using the Galerkin method, in which the mass normalized mode shapes of the simply supported uniform bare plate are used as the basis functions for the out-of-plane displacement. The term “bare plate” here refers to the plate without the random point masses and springs. The first 927 uniform bare plate modes are retained for the analysis which means that all the uniform bare plate modes below 5000 Hz are included. This guarantees the convergence of the energy density at the maximal frequency of interest for all the simulated cases (less than 2 percent difference in the energy density response when the number of uniform plate modes is increased by 20 percent, i.e. from 927 to 1117). For the construction of the system matrices, the bare plate modes are discretized using 263 points along the x1-axis and 315 points along the x2-axis. This corresponds to a mesh density of 15 points per bending wavelength of the uniform plate at 5000 Hz.

3.2. Mean energy response under rain-on-the-roof forcing

3 2 1 0 0.2 0.4

x (m ) 1

0.6 0.8

0.2

0.4

0.6

0.8

m) x 2(

1

Normalized mean squared displacement

Normalized mean squared displacement

Figs. 3–5 present the empirical results of the plates obtained for the mean of the squared displacement jWðxÞj2 , of the ~ energy density eðxÞ and of the local modal energy eðxÞ=nðxÞ. In all cases rain-on-the-roof forcing has been applied and the frequency of interest is 1000 Hz. The squared displacement, energy density and local modal energy have been computed using Eqs. (35), (5) and (11). To present the results, 51 points along the x1-axis and 61 points along the x2-axis have been used to discretize the plates. However, to avoid showing response concentration factors near the edges of the plates [1, Section 13.3], points ðx1 ; x2 Þ outside the range ð½0:1L1 ; 0:9L1 ; ½0:1L2 ; 0:9L2 Þ are not shown in the figures. Moreover, the results are normalized by their spatially averaged value and are smoothed by averaging the value at each point with the values at the eight neighboring points. This light smoothing makes the general trends clearer without averaging them out. As expected Figs. 3 and 4 show that the absolute squared displacement and the energy density of the plates are higher where the plate is thinner while Fig. 5 shows that the local modal energy is spatially uniform. The average deviation of the squared displacement from its spatial mean is 49 percent for the tapered plate and 37 percent for the bumped plate while the average deviation of the energy density is 26 percent for the tapered plate and 22 percent for the bumped plate. The fact that the energy density deviates less than the squared displacement is consistent with Eqs. (5), (10) and (11) since if all the 2 1 properties of the plate are uniform except the thickness hðxÞ, then E½jWðxÞj2  p h ðxÞ while E½eðxÞ p h ðxÞ. In Fig. 5a small but systematic oscillations in the local modal energy are observed near the thicker edge of the tapered plate. These oscillations are not a convergence artefact; they are caused by the fact that certain modes of the system localize in the thinner region of the plate. This feature is illustrated in Fig. 6 where nine mode shapes of the bare tapered plate with natural frequency fr close to 1000 Hz are plotted. It can be seen that modes with a high ratio k2 =k1 localize in the thinner region of the plate (total internal reflection occurs). These are the so-called “wedge modes” of the plate [28]. Because these modes do not propagate in the thicker region of the plate, the local modal energy lacks sufficient wave components to be fully flat in that region. This causes the small oscillations observed in Fig. 5a. In Fig. 5b similar oscillations are observed near the center of the bumped plate. Nevertheless, the oscillations are small for both plates (on average 5–10 percent), and thus Eq. (22) remains applicable.

3 2 1 0 0.2 0.4 x (m ) 1

0.6 0.8

0.2

0.4

0.6

0.8

1

m) x 2(

Fig. 3. Squared displacement response at 1000 Hz under rain-on-the-roof forcing: (a) Tapered plate; (b) bumped plate.

7226

J. Legault et al. / Journal of Sound and Vibration 333 (2014) 7216–7232

3 Normalized mean energy density

Normalized mean energy density

3

2

1

2

1

0

0 0.2

0.4

x (m ) 1

0.6

0.8

0.2

0.4

0.2

1

0.8

0.6

0.4 x ( 1 m)

)

x 2 (m

0.6 0.8

0.2

0.4

1

0.8

0.6

) x 2 (m

3

Normalized mean local modal energy

Normalized mean local modal energy

Fig. 4. Energy density response at 1000 Hz under rain-on-the-roof forcing: (a) Tapered plate; (b) bumped plate.

2

1

0 0.2

0.4 0.6 x (m 1 )

0.8

0.2

0.4

1

0.8

0.6

3

2

1

0 0.2 0.4 0.6 x ( 1 m)

) x 2 (m

0.8

0.2

0.4

0.6

1

0.8

x 2 (m)

Amp. (kg−1/2)

0.5

1

x 1 (m)

2 0 −2 1

0.5

0 0

0.5

1

x 1 (m)

2 0 −2 1

0.5

0 0

0.5

x 1 (m)

0.5

1

0 0

0.5

1

x 1 (m)

2 0 −2 1

0.5

0 0

0.5

1

x 1 (m)

2 0 −2 1

0.5

0 0

0.5

x 1 (m)

2 0 −2 1

0.5

1

0 0

0.5

1

x 1 (m)

fr = 995.09 Hz 2 0 −2 1

0.5

x (m 2 )

fr = 1004.48 Hz

x ( 2 m)

fr = 979.34 Hz

x (m 2 )

fr = 990.00 Hz

x (m ) 2

fr = 1001.19 Hz

x (m 2 )

1

x (m ) 2

fr = 986.70 Hz

x (m ) 2 Amp. (kg−1/2)

0 0

0 −5

Amp. (kg−1/2)

0.5

Amp. (kg−1/2)

1

x ( 2 m)

5

Amp. (kg−1/2)

0 −2

fr = 979.17 Hz

Amp. (kg−1/2)

2

Amp. (kg−1/2)

fr = 978.80 Hz

Amp. (kg−1/2)

Amp. (kg−1/2)

Fig. 5. Local modal energy response at 1000 Hz under rain-on-the-roof forcing: (a) Tapered plate; (b) bumped plate.

0 0

0.5

1

x 1 (m)

fr = 1004.61 Hz 2 0 −2 1

0.5

x (m 2 )

0 0

0.5

x 1 (m)

Fig. 6. Amplitude of nine bare tapered plate mode shapes with natural frequency in the vicinity of 1000 Hz.

1

J. Legault et al. / Journal of Sound and Vibration 333 (2014) 7216–7232

7227

To present a situation where the breakdown of Eq. (22) is more significant, the local modal energy of a “beamed” plate, i.e. a plate with a one-dimensional bump, has been computed. The thickness profile of the plate is shown in Fig. 7a and the local modal energy results are shown in Fig. 7b. Similar to the bumped plate, a Gaussian profile has been used for the onedimensional bump, and the mass and peak ratios of the beam were set to 0.2 and 3. In Fig. 7b, systematic deviations of up to 20 percent are observed for the normalized local modal energy. Whether the plate should still be treated as a single subsystem is thus becoming questionable. If the peak ratio of the one-dimensional bump was increased further, the system would gradually tend towards a built-up structure comprising two plates coupled through a beam. In that case, splitting the structure into three separate subsystems would be necessary. These observations raise the following question: could a significant breakdown of Eq. (22) be used to indicate when a system needs to be split into separate subsystems? This suggestion would merit attention in future investigations. Before ending this section, it should be noted that although plates with spatially varying thickness are considered in this study, Eq. (22) has also been validated for plates with spatially varying Young's modulus or mass per unit volume. The results are not presented here to avoid redundancy. 3.3. Mode shape statistics In Fig. 8, the empirical mean modal norm E½ϕr ðxÞ of the tapered plate and the bumped plate is compared with the predicted value given by Eq. (27). Two locations are considered for the tapered plate: one in the thinner region of the plate, point ð0:23L1 ; 0:21L2 ), and one in the thicker region, point ð0:79L1 ; 0:77L2 Þ. For the bumped plate point ð0:50L1 ; 0:50L2 Þ that corresponds to the center of the bump is considered instead of point ð0:79L1 ; 0:77L2 Þ. In the figure, the predicted value for the uniform plate ð ¼ 1=M 0 Þ is also shown for comparison. To compute the mean modal norm at frequency ω, the mode shape ϕr ðxÞ with natural frequency ωr closest to ω has been found for each realization of the system, and the ensemble average has been performed over the obtained 5000 samples. Good agreement is observed between the empirical and the predicted results. Using the same 5000 samples, the empirical mode shape factor KðxÞ has also been computed and is plotted in Fig. 9a. For both plates, KðxÞ tends towards 3 with increasing frequency, which suggests that the mode shapes are becoming Gaussian, and hence that GOE statistics also apply to the mode shapes of inhomogeneous systems with slowly varying properties. To explore this suggestion, the pdf of the mode shape with natural frequency closest to 1000 Hz is shown in Fig. 9b. The mode shape amplitude has been normalized by its root mean square value, so that the pdf relates to the variable

Normalized mean local modal energy

Normalized thickness

2

3 2 1 0 0

2 1 0 0.2

1 0.5 x ( 1 m)

3

0.4

0.5 1

x ( 1 m)

) x 2 (m

0

0.6 0.8

0.4

0.2

0.6

1

0.8

) x 2 (m

Fig. 7. Beamed plate plots: (a) Normalized thickness; (b) local modal energy at 1000 Hz under rain-on-the-roof forcing.

0.6

0.6

Empirical (0.23L ,0.21L )

Empirical (0.23L ,0.21L )

Predicted (0.23L ,0.21L ) Predicted (0.79L ,0.77L )

0.4

Predicted (uniform plate)

0.3 0.2 0.1 0

Empirical (0.50L ,0.50L ) −1

0.5

Mean modal norm (kg )

−1

Mean modal norm (kg )

Empirical (0.79L ,0.77L )

0.5

Predicted (0.23L ,0.21L ) Predicted (0.50L ,0.50L )

0.4

Predicted (uniform plate)

0.3 0.2 0.1

0

200

400

600

Frequency (Hz)

800

1000

0

0

200

400

600

Frequency (Hz)

Fig. 8. Mean modal norm: (a) Tapered plate; (b) bumped plate.

800

1000

7228

J. Legault et al. / Journal of Sound and Vibration 333 (2014) 7216–7232

0.8

5

Gaussian pdf Tapered plate (0.23L1,0.21L2)

0.7

3

2 Tapered plate (0.23L1,0.21L2)

Probability density

Mode shape factor

4

Bumped plate (0.23L1,0.21L2)

0

200

400

600

800

Bumped plate (0.50L1,0.50L2)

0.4 0.3

0.1

Bumped plate (0.50L1,0.50L2)

0

Bumped plate (0.23L1,0.21L2)

0.5

0.2

Tapered plate (0.79L1,0.77L2)

1

Tapered plate (0.79L1,0.77L2)

0.6

0 −3

1000

−2

Frequency (Hz)

−1

0

1

2

3

Normalized mode shape

Mean of real part of input mobility (ms

Empirical (0.23L ,0.21L ) Empirical (0.79L ,0.77L ) Predicted (0.23L ,0.21L )

10

Predicted (0.79L ,0.77L ) Predicted (uniform plate)

10

10

10

0

200

400

600

Frequency (Hz)

800

1000

N

N

10

Mean of real part of input mobility (ms

)

)

Fig. 9. Higher order statistics of the mode shapes: (a) Mode shape factor; (b) normalized mode shape pdf at 1000 Hz.

10

Empirical (0.23L ,0.21L ) Empirical (0.50L ,0.50L ) Predicted (0.23L ,0.21L )

10

Predicted (0.50L ,0.50L ) Predicted (uniform plate)

10

10

10

0

200

400

600

800

1000

Frequency (Hz)

Fig. 10. Ensemble average of the real part of the input mobility: (a) Tapered plate; (b) bumped plate.

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

ϕr ðxÞ= E½ϕ2r ðxÞ. In the figure, the empirical pdfs are compared with the pdf of a standard Gaussian variable. Excellent agreement is obtained between the empirical pdfs and the Gaussian pdf at all locations.

3.4. Mean input mobility Empirical results for the ensemble average of the real part of the input mobility are compared in Fig. 10 with the predicted value given by Eq. (32). The locations considered in the previous subsection for the modal norm are also considered here. In the figure, the predicted value for the uniform plate is also shown for comparison. Good agreement is obtained between the empirical results and the predictions of Eq. (32). In contrast, the uniform plate prediction overestimates the mobility in the thicker region of the tapered plate by a factor of 2 while it underestimates the mobility in the thinner region by slightly less than a factor of 2. For the bumped plate, it overestimates the mobility at the center of the bump by a factor of 10. These differences justify the use of Eq. (32) and underline the pertinence of this study. The good agreement obtained in Fig. 10 also validates the formula for the ensemble-averaged total energy under point forcing since there is only a deterministic factor jF 0 j2 =2ηω between the real part of the input mobility and the total energy (see Section 2.4).

3.5. Relative variance of the total energy Results for the ensemble relative variance of the total energy of the tapered plate under rain-on-the-roof and point forcing are shown in Fig. 11. These results are compared with the predictions of Eq. (29) using α ¼1 for rain-on-the-roof forcing and α ¼ KðxÞ ¼ 3 for point forcing. Good agreement is obtained in both cases. Good agreement has also been obtained for the bumped plate but to avoid redundancy the results are not reproduced here. For rain-on-the-roof forcing the relative variance of the energy is only due to the randomness of the natural frequencies. Therefore the good agreement obtained in Fig. 11 confirms that GOE statistics also apply to the statistics of the natural frequencies of inhomogeneous systems with slowly varying properties. For point forcing, the relative variance is slightly underestimated by the predicted value at point ð0:23L1 ; 0:21L2 Þ but it is well estimated at point ð0:79L1 ; 0:77L2 Þ. This is because

J. Legault et al. / Journal of Sound and Vibration 333 (2014) 7216–7232

Relative variance of total energy

10

7229

Empirical: point forcing, (0.23L ,0.21L ) Empirical: point forcing, (0.79L ,0.77L ) Empirical: rain−on−the−roof forcing

10

Predicted: point forcing Predicted: rain−on−the−roof forcing

10

10

10

0

200

400

600

800

1000

Frequency (Hz)

Fig. 11. Ensemble relative variance of the total energy of the tapered plate.

1

0.5

0

−0.5 Empirical Predicted: corrected phase Predicted: uncorrected phase −1 0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

Real part of normalized TF

Imaginary part of normalized TF

1

0.5

0

−0.5

−1 0.1

Empirical Predicted: corrected phase Predicted: uncorrected phase 0.2

0.3

0.4

x1 (m)

0.5

0.6

0.7

0.8

x1 (m)

Fig. 12. Normalized transfer function of the tapered plate at 1000 Hz with x0 ¼ ð0:50L1 ; 0:50L2 Þ: (a) Imaginary part; (b) real part.

the modal coefficient KðxÞ is close to 3 at point ð0:79L1 ; 0:77L2 Þ while it is slightly larger at point ð0:23L1 ; 0:21L2 Þ (see Fig. 9a). Such spatial variations of the mode shape factor within a system have also been reported by Langley and Brown [14]. 3.6. Transfer function statistics To assess the accuracy of the corrected phase conjecture made in Section 2.5, the normalized transfer function, defined as Hnorm ðx; x0 Þ ¼

 E½Hðx; x0 Þ ðE½ImfHðx; xÞgE½ImfHðx0 ; x0 ÞgÞ1=2

;

(55)

has been computed. If the phase conjecture is correct, it follows from Eqs. (39) and (47) that the imaginary and real parts of the normalized transfer function will be given by ImfHnorm ðx; x0 Þg ¼ RefΨ ðfkRgcorr Þg ¼ J 0 ðfkRgcorr Þ;

(56a)

RefH norm ðx; x0 Þg ¼ ImfΨ ðfkRgcorr Þg:

(56b)

It should be noted that the imaginary part of H norm ðx; x0 Þ corresponds to the diffuse field correlation function Cðx; x0 Þ. The empirical results for normalized transfer function are shown in Fig. 12 for the tapered plate and in Fig. 13 for the bumped plate. The center point of the plate has been chosen as the drive location x0 and 201 observation locations x have been considered along the x2 ¼ 0:50L2 line. The frequency of interest is 1000 Hz. In Figs. 12 and 13, the empirical normalized transfer function is compared with two predictions: the corrected and the uncorrected phase prediction. The corrected phase prediction corresponds to the prediction of Eq. (56) while the uncorrected phase prediction is calculated by replacing the corrected phase with the uncorrected phase in Eq. (56). The uncorrected phase is obtained by multiplying the wavenumber k at the reference location with the distance R between the reference and the observation location. Comparing the corrected and uncorrected phase predictions allows us to assess if using the corrected phase makes a worthwhile difference. For the tapered plate, relatively good agreement is obtained between the empirical result and the corrected phase prediction (the asymmetry of the plate is well captured). In comparison, the uncorrected phase prediction gives less good

7230

J. Legault et al. / Journal of Sound and Vibration 333 (2014) 7216–7232

1

0.5

0

−0.5

−1 0.1

Empirical Predicted: corrected phase Predicted: uncorrected phase 0.2

0.3

0.4

0.5

0.6

0.7

Real part of normalized TF

Imaginary part of normalized TF

1

0.5

0

−0.5

−1 0.1

0.8

Empirical Predicted: corrected phase Predicted: uncorrected phase 0.2

0.3

0.4

x1 (m)

0.5

0.6

0.7

0.8

x1 (m)

Fig. 13. Normalized transfer function of the bumped plate at 1000 Hz with x0 ¼ ð0:50L1 ; 0:50L2 Þ: (a) Imaginary part; (b) real part.

1

0.5

0

−0.5 Empirical Predicted: corrected phase Predicted: uncorrected phase −1 0.1

0.2

0.3

0.4

0.5 x (m)

0.6

0.7

0.8

Real part of normalized TF

Imaginary part of normalized TF

1

0.5

0

−0.5

−1 0.1

Empirical Predicted: corrected phase Predicted: uncorrected phase 0.2

0.3

0.4

0.5

0.6

0.7

0.8

x (m)

Fig. 14. Normalized transfer function of the bumped plate at 1000 Hz with x0 ¼ ð0:35L1 ; 0:50L2 Þ: (a) Imaginary part; (b) real part.

results. For x1 40:80 m, the amplitude of the imaginary does not have the form of the Bessel function J0, and therefore neither the corrected phase prediction nor the uncorrected phase prediction yields good agreement in that range. This discrepancy may be related to the lack of wave components mentioned in Section 3.2 although more investigations are needed to clarify this issue. For the bumped plate, relatively good agreement is obtained for the phase of the normalized transfer function (the empirical and predicted peaks and dips are approximately aligned), but the amplitude of the real and imaginary parts is less well captured. The fact that the phase is approximately captured is not surprising since the center point of the plate also corresponds to the center of the axisymmetric bump, and therefore the corrected prediction ought to work approximately, at least for the phase. A case that is likely to generate more distortion is when the drive point x0 is offset from the center of the bump. Results for such a case are presented in Fig. 14. The drive point x0 is ð0:35L1 ; 0:50L2 Þ. For x1 4 0:55 m, the empirical normalized transfer function does not have the form of the cylindrical wave function Ψ , and therefore the predictions fail in that range. This illustrates the limits of the conjecture made in Section 2.5. However, it should be remarked that failing to predict the direct field component far from the drive point may not be critical when predicting the mean squared transfer function magnitude. Indeed, it has been shown in Section 2.5 that far from the drive point, the reverberant field will dominate the response (see Eq. (50)), and hence errors in the prediction of the direct field will not make an important difference. To illustrate this, Fig. 15 presents the empirical results obtained for E½jHðx; x0 Þj2 . The empirical results are compared with the prediction obtained from Eqs. (40), (47) and (49). For the tapered plate, E½jHðx; x0 Þj2  has been computed at points x ¼ ð0:55L1 ; 0:50L2 Þ and x ¼ ð0:85L1 ; 0:50L2 Þ, and the point force has been applied at point x0 ¼ ð0:50L1 ; 0:50L2 Þ. The average trend of E½jHðx; x0 Þj2  is relatively well captured by the predictions at both points. Because point ð0:55; 0:50L2 Þ is close to the drive point, a large error in the prediction of the direct field would have yielded a large error between the predicted and the empirical squared displacement. On the other hand, point ð0:85L1 ; 0:50L2 Þ is far from the drive point. Therefore, even if the direct field component is not well predicted at that point (see Fig. 12), the predicted and empirical results are in good agreement. For the bumped plate, E½jHðx; x0 Þj2  has been computed at points x ¼ ð0:40L1 ; 0:50L2 Þ and x ¼ ð0:80L1 ; 0:50L2 Þ, and the point force has been applied at point x0 ¼ ð0:35L1 ; 0:50L2 Þ. Good agreement is obtained at both locations, which corroborates the above analysis.

J. Legault et al. / Journal of Sound and Vibration 333 (2014) 7216–7232

10 Empirical (0.55L ,0.50L ) Empirical (0.85L ,0.50L ) Predicted (0.55L ,0.50L ) Predicted (0.85L ,0.50L )

10

10

0

200

400

600

Frequency (Hz)

800

1000

Mean squared TF magnitude (m /N )

Mean squared TF magnitude (m /N )

10

10

7231

Empirical (0.40L ,0.50L ) Empirical (0.80L ,0.50L ) Predicted (0.40L ,0.50L ) Predicted (0.80L ,0.50L )

10

10

10

0

200

400

600

800

1000

Frequency (Hz)

Fig. 15. Mean squared transfer function magnitude: (a) Tapered plate, x0 ¼ ð0:50L1 ; 0:50L2 Þ; (b) bumped plate, x0 ¼ ð0:35L1 ; 0:50L2 Þ.

4. Conclusion In this paper closed-form expressions have been derived, within the spirit of classical statistical energy analysis, for the spatial variation of energy within a single inhomogeneous system with slowly varying properties. The key result has been derived by a modal argument based on the strong-coupling limit of SEA, and also by a wave argument based on energy flux conservation in the absence of significant reflections within the system. The main results are Eq. (13) for the scaling of the mean energy density, but other SEA-like results have also been obtained: Eq. (27) for the mean modal norm and Eq. (32) for the mean value of the real part of the input mobility. As mentioned in the Introduction, other approaches could be used to predict the scaling of the energy density within an inhomogeneous system, but these approaches are more complex and they would yield numerical predictions, not closedform expressions. Hence this study represents a step forward to obtain local information inside an inhomogeneous SEA subsystem. Numerical studies have been presented, which suggest that the analysis holds to a reasonable degree of accuracy even when the spatial inhomogeneity is surprisingly large. More speculative but also useful is Eq. (46), for the corrected phase of the direct field Green's function. This approximate treatment has been shown to work reasonably well in some circumstances, but more investigations are required to determine the conditions under which it is valid. It has been assumed in this work that the system has a single scalar response variable, such as the out-of-plane displacement for a plate. The extension of the theory to inhomogeneous systems with more than one response variable, such as plates in which the in-plane motion is coupled to the out-of-plane motion due to spatially-varying curvature, seems straightforward in principle but has yet to be tested in detail. Moreover, the present analysis has not considered the situation when the inhomogeneous subsystem is coupled to other subsystems. Developing a fully-fledged version of SEA for inhomogeneous subsystems will require the calculation of coupling loss factors between inhomogeneous subsystems. It is easy to see how such a calculation could be approached using similar approximations and assumptions to the theory presented here, but the implementation and testing of such a calculation is a topic for future work. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17]

R.H. Lyon, R.G. DeJong, Theory and Application of Statistical Energy Analysis, 2nd edition, Butterworth-Heinemann, Newton, MA, 1995. B.R. Mace, Statistical energy analysis, energy distribution models and system modes, Journal of Sound and Vibration 264 (2003) 391–409. B.R. Mace, Statistical energy analysis: coupling loss factors, indirect coupling and system modes, Journal of Sound and Vibration 279 (2005) 141–170. A.D. Pierce, Physical interpretation of the wkb or Eikonal approximation for waves and vibrations in inhomogeneous beams and plates, The Journal of the Acoustical Society of America 48 (1970) 275–284. O.A. Germogenova, Geometrical theory for flexure waves in shells, The Journal of the Acoustical Society of America 53 (1973) 535–540. M.J. Lighthill, Energy flow in the cochlea, Journal of Fluid Mechanics 106 (1981) 149–213. D.J. O'Boy, V. Krylov, V. Kralovic, Damping of flexural vibrations in rectangular using the acoustic black hole effect, Journal of Sound and Vibration 329 (2010) 4672–4688. E. Savin, High-frequency dynamics of heterogeneous slender structures, Journal of Sound and Vibration 332 (2013) 2461–2487. N. Totaro, J.L. Guyader, Extension of the statistical modal energy distribution analysis for estimating energy density in coupled subsystems, Journal of Sound and Vibration 331 (2012) 3114–3129. G. Tanner, Dynamical energy analysis—determining wave energy distributions in vibro-acoustical structures in the high-frequency regime, Journal of Sound and Vibration 320 (2009) 1023–1038. F.J. Fahy, P. Gardonio, Sound and Structural Vibration, 2nd edition, Elsevier, Oxford, UK, 2007. M.L. Mehta, Random Matrices, 2nd edition, Academic Press, Boston, USA, 1991. R. Weaver, Spectral statistics in elastodynamics, The Journal of the Acoustical Society of America 85 (1989) 1005–1013. R.S. Langley, A.W.M. Brown, The ensemble statistics of the energy of a random system subjected to harmonic excitation, Journal of Sound and Vibration 275 (2004) 823–846. P.W. Smith, Statistical models of coupled dynamical systems and the transition from weak to strong coupling, The Journal of the Acoustical Society of America 65 (1979) 695–698. M.V. Berry, Semi classical mechanics of regular and irregular motion, in: G. Jooss, R.H.G. Helleman, R. Stora (Eds.), Les Houches, Lecture Series Session, vol. XXXVI, 1983, pp. 171–271. R.S. Langley, The modal density of anisotropic structural components, The Journal of the Acoustical Society of America 99 (1996) 3481–3487.

7232

J. Legault et al. / Journal of Sound and Vibration 333 (2014) 7216–7232

[18] R.S. Langley, A wave intensity technique for the analysis of high frequency vibrations, Journal of Sound and Vibration 159 (1992) 483–502. [19] R.S. Langley, Elastic wave transmission coefficients and coupling loss factors for structural junctions between curved panels, Journal of Sound and Vibration 169 (1994) 297–317. [20] O.A. Germogenova, Formalism of geometrical optics for flexure waves, The Journal of the Acoustical Society of America 49 (1971) 776–780. [21] R.H. Lyon, Needed: a new definition of diffusion, The Journal of the Acoustical Society of America 56 (1974) 1300–1302. [22] R.S. Langley, V. Cotoni, The ensemble statistics of the vibrational energy density of a random system subjected to single point harmonic excitation, The Journal of the Acoustical Society of America 118 (2005) 3064–3076. [23] L. Cremer, M. Heckl, B.A. Petersson, Structure-Borne Sound: Structural Vibrations and Sound Radiation at Audio Frequencies, 3rd edition, ButterworthHeinemann, Newton, MA, 2005. [24] O.I. Lobkis, R.L. Weaver, On the emergence of the green's function in the correlations of a diffuse field, The Journal of the Acoustical Society of America 110 (2001) 3011–3017. [25] P.J. Shorter, R.S. Langley, On the reciprocity relationship between direct field radiation and diffuse reverberant loading, The Journal of the Acoustical Society of America 117 (2005) 85–95. [26] W. Choi, R.S. Langley, J. Woodhouse, Boundary effects on the vibration statistics of a random plate, Journal of Sound and Vibration 332 (2013) 850–866. [27] F.C. Appl, N.R. Byers, Fundamental frequency of simply supported rectangular plates with linearly varying thickness, Journal of Applied Mechanics 32 (1965) 163–168. [28] V.V. Krylov, Geometrical-acoustics approach to the description of localized vibrational modes of an elastic solid wedge, Soviet Physics: Technical Physics 35 (1990) 137–140.