Vibrational energy flow models for dilatational wave in elastic solids

Abstract

The analysis of medium- to high-frequency vibrations of structures is of particular interest in various fields, including the aviation and ship-building industries. Energy flow analysis is well known to be effective for structural acoustics problems in the medium- to high-frequency range, with sufficient detail to include all significant information of the vibrational energy levels or global variation in vibrational energy densities. The aim of this article is to develop energy flow models for a dilatational wave in three-dimensional elastic solids. The energy governing equation derived for the model is expressed in terms of the time- and space-averaged energy density, which represents the global variation of the energy density quite well. Numerical analyses are performed to verify the validity and accuracy of the model for a cube-shaped structure vibrating at a single frequency, and the results of the analysis of energy density distributions from the energy flow analysis that is developed are compared to those obtained using NASTRAN. The energy flow models for the elastic solids, which are useful in the prediction of the vibrational response for a three-dimensional structural analysis at a medium-to-high frequency, are newly derived.

Introduction

During the design phase, a vibrational analysis is implemented for built-up structures to predict the vibrational level of the structures and to identify structure-borne noise paths. Over the last few decades, advances in technology have drawn attention to high-frequency noise and vibration. Although finite element analysis (FEA) or boundary element analysis (BEA) is often used to predict vibrational response, FEA and BEA are difficult to apply for high-frequency vibrational analyses because the mesh size requirement is unreasonably small, resulting in a high computational intensity.

Various approaches have been attempted in order to overcome the limitations of conventional methods. For example, statistical energy analysis (SEA) is a representative method using statistical approaches, and it has been used to predict the space- and frequency-averaged response of the structure at high frequencies. However, SEA has some limitations in offering detailed information, such as the spatial variation of energy density or energy flow that are essential in the design of structures. As an alternative, energy flow analysis (EFA) was developed in order to predict the vibrational response of the structure. EFA formulates the flow of mechanical energy, which is analogous to heat conduction. Since the energy governing equation for EFA is in the form of a parabolic partial differential equation, the solution can suggest the spatial variation of energy densities or energy flow paths.

The EFA was first introduced by Belov et al. ¹ Nefske and Sung ² applied FEA to the energy governing equation in order to predict the vibrational response of the Euler–Bernoulli beam. Wohlever and Bernhard ³ established the energy flow model based on classical displacement solutions for a harmonically excited rod and Euler–Bernoulli beam. Bouthier and Bernhard^4,5 also derived energy governing equations for the membrane and Kirchhoff plate. Cho ⁶ developed the joint relationships for coupled structures based on the concept of a power transmission and reflection coefficient. Park et al.^7,8 developed energy flow models for the in-plane waves of isotropic thin plates, as well as a flexural wave in an orthotropic thin plate. Park and Hong^9–11 developed the energy governing equations of the Timoshenko beam and Mindlin plate. Song and Hong ¹² developed a non-conservative joint for the EFA of a coupled beam structure, and Song et al. ¹³ extended this model for penetration of beam-plate structures. Kwon et al. ¹⁴ derived the energy governing equation for the out-of-plane wave in a doubly curved shell with two curvatures along each orthogonal coordinate axes (Table 1).

Table 1.

Summary of related studies and comparison with this study.

Studies	Year	Target structure	Wave type
Nefske and Sung 2	1989	Euler–Bernoulli beam	Flexural wave
Bouthier and Bernhard4,5	1995	Membrane	Out-of-plane wave
		Isotropic thin plate	Flexural wave
Park et al.7,8	2001	Isotropic thin plate	In-plane wave
	2003	Orthotropic thin plate	In-plane wave
Park and Hong9–11	2006	Timoshenko beam	Flexural wave
	2008	Mindlin plate	Out-of-plane wave
Kwon et al. 14	2012	Doubly curved shell	Out-of-plane wave
This study	2016	Isotropic solid	Dilatational wave

In prior EFA, the energy governing equation was derived for a one-dimensional (1D) or two-dimensional (2D) structure, such as a beam, plate, membrane, and shell. However, there are many industrial situations in vibration analysis where it is necessary to determine the three-dimensional (3D) vibrational response, such as in the pier of a bridge, core of an industrial generator, or in an automotive frame. Usually, these structures can be analyzed by modeling with a solid element based on finite element method (FEM) at low frequencies. Even then, it is impossible to apply these to the vibrational problem at a medium-to-high frequency. Since the precise vibrational analysis for 3D structure vibrating at a medium-to-high frequency is required, the development of an energy flow model for a 3D structure is highly desirable.

In 3D elastic solids, structural waves can be composed of compressional waves and shear waves, which are referred to as dilatational waves and distortion waves, respectively. The aim of this work is to develop an energy flow model for a dilatational wave in 3D elastic solids vibrating at medium- to high-frequency ranges. Numerical analyses are implemented for a cube-shaped solid excited by a spherical cavity source, which is subjected to the internal pressure. To verify the derived model, the results of the analysis obtained from the derived EFA model are compared to those from NASTRAN, and the effects of the structural damping and frequency are investigated.

Derivation of the energy governing equations for a dilatational wave in an elastic solid medium

Unlike beam or plate structures, there is no obvious approximate theory for the kinematics of deformation. Thus, the energy governing equation for a solid element can be derived from the exact elasticity equation. The exact governing equation for a homogeneous isotropic elastic solid in terms of displacements is ¹⁵

$$\mu {\nabla }^2\overrightarrow{\mathrm{u}}+(\lambda +\mu )\nabla \nabla ·\overrightarrow{\mathrm{u}}=\rho \frac{{\partial }^2\overrightarrow{\mathrm{u}}}{\partial {\mathrm{t}}^2}$$

(1)

where $\overrightarrow{\mathrm{u}}$ is the particle displacement vector, $\rho$ is the mass density of solid, $\lambda$ is Lame’s first constant, and $\mu$ is Lame’s second constant. The displacement vector can be composed by two potential functions using the Helmholtz resolution theorem as ¹⁶

$$\overrightarrow{\mathrm{u}}=\nabla \phi +\nabla \times \overrightarrow{\psi }$$

(2)

where $\phi$ is a scalar potential function that describes the change in volume or dilatational disturbance. $\overrightarrow{\psi }$ is a vector potential function that represents the rotation or distortional disturbance of elastic solid. Substituting equation (2) into equation (1) yields

$$\nabla [(\lambda +2\mu ){\nabla }^2\phi -\rho \frac{{\partial }^2\phi }{\partial {\mathrm{t}}^2}]+\nabla \times [\mu {\nabla }^2\overrightarrow{\psi }-\rho \frac{{\partial }^2\overrightarrow{\psi }}{\partial {\mathrm{t}}^2}]=0$$

(3)

Equation (3) will be satisfied if each bracketed term vanishes. The first bracket yields the governing equation for the dilatational wave in terms of the scalar potential

$${\nabla }^2\phi =\frac{1}{{\mathrm{c}}_{\mathrm{L}}^2}\frac{{\partial }^2\phi }{\partial {\mathrm{t}}^2}$$

(4)

$${\mathrm{c}}_{\mathrm{L}}={\left[\frac{\lambda +2\mu }{\rho }\right]}^{1/2}$$

(5)

where ${\mathrm{c}}_{\mathrm{L}}$ is the wave speed of the dilatational wave. The second bracket yields the governing equation for the distortion wave in terms of the vector potential

$${\nabla }^2\overrightarrow{\psi }=\frac{1}{{{\mathrm{c}}_{\mathrm{T}}}^2}\frac{{\partial }^2\overrightarrow{\psi }}{\partial {\mathrm{t}}^2}$$

(6)

$${\mathrm{c}}_{\mathrm{T}}={\left[\frac{\mu }{\rho }\right]}^{1/2}$$

(7)

where ${\mathrm{c}}_{\mathrm{T}}$ is the wave speed of the distortion wave.

Since this article covers only the dilatational wave, it is assumed that there is no distortion wave type, and the displacement in the elastic solid can thus be described only with the scalar potential function. In the elastic solid without a distortion wave, the time-averaged potential energy density $〈{\mathrm{e}}_{\mathrm{p}}〉$ , time-averaged kinetic energy density $〈{\mathrm{e}}_{\mathrm{k}}〉$ , and time-averaged intensity $〈{\mathrm{I}}_{\mathrm{i}}〉$ can be derived in terms of the scalar potential function based on the stress–strain relationship as

$$〈{\mathrm{e}}_{\mathrm{p}}〉=\frac{1}{2}\mathrm{Re}[+\mu \begin{array}{c}\frac{\lambda }{2}{\nabla }^2\phi ({\nabla }^2\phi {)}^{*}\hfill \\ \multicolumn{1}{c}{\left\{\begin{array}{c}\frac{{\partial }^2\phi }{\partial {\mathrm{x}}^2}{\left(\frac{{\partial }^2\phi }{\partial {\mathrm{x}}^2}\right)}^{*}+\frac{{\partial }^2\phi }{\partial {\mathrm{y}}^2}{\left(\frac{{\partial }^2\phi }{\partial {\mathrm{y}}^2}\right)}^{*}\hfill \\ +\frac{{\partial }^2\phi }{\partial {\mathrm{z}}^2}{\left(\frac{{\partial }^2\phi }{\partial {\mathrm{z}}^2}\right)}^{*}+\frac{{\partial }^2\phi }{\partial \mathrm{x}\partial \mathrm{y}}{\left(\frac{{\partial }^2\phi }{\partial \mathrm{x}\partial \mathrm{y}}\right)}^{*}\hfill \\ +\frac{{\partial }^2\phi }{\partial \mathrm{y}\partial \mathrm{z}}{\left(\frac{{\partial }^2\phi }{\partial \mathrm{y}\partial \mathrm{z}}\right)}^{*}+\frac{{\partial }^2\phi }{\partial \mathrm{z}\partial \mathrm{x}}{\left(\frac{{\partial }^2\phi }{\partial \mathrm{z}\partial \mathrm{x}}\right)}^{*}\hfill \end{array}\right\}}\end{array}]$$

(8)

$$〈{\mathrm{e}}_{\mathrm{k}}〉=\frac{1}{2}\mathrm{Re}\left[\frac{1}{2}\rho \left\{\begin{array}{c}\frac{{\partial }^2\phi }{\partial \mathrm{x}\partial \mathrm{t}}{\left(\frac{{\partial }^2\phi }{\partial \mathrm{x}\partial \mathrm{t}}\right)}^{*}+\frac{{\partial }^2\phi }{\partial \mathrm{y}\partial \mathrm{t}}{\left(\frac{{\partial }^2\phi }{\partial \mathrm{y}\partial \mathrm{t}}\right)}^{*}\hfill \\ \multicolumn{1}{c}{+\frac{{\partial }^2\phi }{\partial \mathrm{z}\partial \mathrm{t}}{\left(\frac{{\partial }^2\phi }{\partial \mathrm{z}\partial \mathrm{t}}\right)}^{*}}\end{array}\right\}\right]$$

(9)

$$〈{\mathrm{I}}_{\mathrm{x}}〉=-\frac{1}{2}\mathrm{Re}\left[\begin{array}{c}(\lambda {\nabla }^2\phi +2\mu \frac{{\partial }^2\phi }{\partial {\mathrm{x}}^2}){\left(\frac{{\partial }^2\phi }{\partial \mathrm{x}\partial \mathrm{t}}\right)}^{*}+\hfill \\ \multicolumn{1}{c}{2\mu \{\frac{{\partial }^2\phi }{\partial \mathrm{x}\partial \mathrm{y}}{\left(\frac{{\partial }^2\phi }{\partial \mathrm{y}\partial \mathrm{t}}\right)}^{*}+\frac{{\partial }^2\phi }{\partial \mathrm{x}\partial \mathrm{z}}{\left(\frac{{\partial }^2\phi }{\partial \mathrm{z}\partial \mathrm{t}}\right)}^{*}\}}\end{array}\right]$$

(10)

$$〈{\mathrm{I}}_{\mathrm{y}}〉=-\frac{1}{2}\mathrm{Re}\left[\begin{array}{c}(\lambda {\nabla }^2\phi +2\mu \frac{{\partial }^2\phi }{\partial {\mathrm{y}}^2}){\left(\frac{{\partial }^2\phi }{\partial \mathrm{y}\partial \mathrm{t}}\right)}^{*}+\hfill \\ \multicolumn{1}{c}{2\mu \{\frac{{\partial }^2\phi }{\partial \mathrm{y}\partial \mathrm{x}}{\left(\frac{{\partial }^2\phi }{\partial \mathrm{x}\partial \mathrm{t}}\right)}^{*}+\frac{{\partial }^2\phi }{\partial \mathrm{y}\partial \mathrm{z}}{\left(\frac{{\partial }^2\phi }{\partial \mathrm{z}\partial \mathrm{t}}\right)}^{*}\}}\end{array}\right]$$

(11)

$$〈{\mathrm{I}}_{\mathrm{z}}〉=-\frac{1}{2}\mathrm{Re}\left[\begin{array}{c}(\lambda {\nabla }^2\phi +2\mu \frac{{\partial }^2\phi }{\partial {\mathrm{z}}^2}){\left(\frac{{\partial }^2\phi }{\partial \mathrm{z}\partial \mathrm{t}}\right)}^{*}+\hfill \\ \multicolumn{1}{c}{2\mu \{\frac{{\partial }^2\phi }{\partial \mathrm{z}\partial \mathrm{x}}{\left(\frac{{\partial }^2\phi }{\partial \mathrm{x}\partial \mathrm{t}}\right)}^{*}+\frac{{\partial }^2\phi }{\partial \mathrm{z}\partial \mathrm{y}}{\left(\frac{{\partial }^2\phi }{\partial \mathrm{y}\partial \mathrm{t}}\right)}^{*}\}}\end{array}\right]$$

(12)

where $\mathrm{Re}[]$ represents the real part, while the $〈〉$ operator and the subscript * indicate the time-averaged quantity and complex conjugate, respectively. With the plane wave assumption, the scalar potential function can be written as shown below, using the wave approach to consider all outgoing wave directions

$$\mathrm{\phi }\left(\mathrm{x},\mathrm{y},\mathrm{z},\mathrm{t}\right)={\mathrm{A}}^{---}+{\mathrm{B}}^{+--}+{\mathrm{C}}^{-+-}+{\mathrm{D}}^{--+}+{\mathrm{E}}^{++-}+{\mathrm{F}}^{+-+}+{\mathrm{G}}^{-++}+{\mathrm{H}}^{+++}$$

(13)

where the unknowns, ${[]}^{\pm \pm \pm }([]=\mathrm{A},\mathrm{B},\mathrm{C},\mathrm{D},\mathrm{E},\mathrm{D},\mathrm{F},\mathrm{G},\mathrm{H})$ , are defined as $[]\times \exp (\pm \mathrm{j}{\mathrm{k}}_{\mathrm{x}}\mathrm{x}\pm \mathrm{j}{\mathrm{k}}_{\mathrm{y}}\mathrm{y}\pm \mathrm{j}{\mathrm{k}}_{\mathrm{z}}\mathrm{z}+\mathrm{j}\mathrm{\omega }\mathrm{t})$ . ${\mathrm{k}}_{\mathrm{x}}$ , ${\mathrm{k}}_{\mathrm{y}}$ , and ${\mathrm{k}}_{\mathrm{z}}$ are the $\mathrm{x}$ , $\mathrm{y}$ , and $\mathrm{z}$ components of the complex wavenumber, ${\mathrm{k}}_{\mathrm{L}}$ , respectively. Substituting equation (13) into equation (4), the dispersion relation of the dilatational wave can be expressed as

$${\mathrm{k}}_{\mathrm{x}}^2+{\mathrm{k}}_{\mathrm{y}}^2+{\mathrm{k}}_{\mathrm{z}}^2=\frac{{\omega }^2}{{\mathrm{c}}_{\mathrm{L}}^2}={\mathrm{k}}_{\mathrm{L}}^2$$

(14)

In the case of small hysteretic damping, the wave numbers are well approximated using a Taylor approximation

$$\begin{array}{c}{\mathrm{k}}_{\mathrm{L}}={\mathrm{k}}_{\mathrm{L}1}-\mathrm{j}{\mathrm{k}}_{\mathrm{L}2}\approx {\mathrm{k}}_{\mathrm{L}1}(1-\mathrm{j}\frac{\eta }{2})\hfill \\ {\mathrm{k}}_{\mathrm{x}}={\mathrm{k}}_{\mathrm{x}1}-\mathrm{j}{\mathrm{k}}_{\mathrm{x}2}\approx {\mathrm{k}}_{\mathrm{x}1}(1-\mathrm{j}\frac{\eta }{2})\hfill \\ {\mathrm{k}}_{\mathrm{y}}={\mathrm{k}}_{\mathrm{y}1}-\mathrm{j}{\mathrm{k}}_{\mathrm{y}2}\approx {\mathrm{k}}_{\mathrm{y}1}(1-\mathrm{j}\frac{\eta }{2})\hfill \\ {\mathrm{k}}_{\mathrm{z}}={\mathrm{k}}_{\mathrm{z}1}-\mathrm{j}{\mathrm{k}}_{\mathrm{z}2}\approx {\mathrm{k}}_{\mathrm{z}1}(1-\mathrm{j}\frac{\eta }{2})\hfill \end{array}$$

(15)

Substituting equation (15) into equation (14), the approximated dispersion relation is presented as

$${\mathrm{k}}_{\mathrm{L}1}^2={\mathrm{k}}_{\mathrm{x}1}^2+{\mathrm{k}}_{\mathrm{y}1}^2+{\mathrm{k}}_{\mathrm{z}1}^2=\frac{{\omega }^2}{{\mathrm{c}}_{\mathrm{L}}^2}$$

(16)

Substituting equation (13) into equations (8) and (9), the expressions for the time-averaged energy densities expanded in terms of the constants A, B, …, H can be written as

$$〈{\mathrm{e}}_{\mathrm{p}}〉=\frac{1}{4}\mathrm{Re}\left[\begin{array}{l}\left\{{\mathrm{\lambda }}_{\mathrm{c}}{\left|{\mathrm{k}}_{\mathrm{L}}^2\right|}^2+2{\mu }_{\mathrm{c}}\left({\left|{\mathrm{k}}_{\mathrm{x}}^2\right|}^2+{\left|{\mathrm{k}}_{\mathrm{y}}^2\right|}^2+{\left|{\mathrm{k}}_{\mathrm{z}}^2\right|}^2\right)\right\}\\ \times {\left|\begin{array}{c}{\mathrm{A}}^{---}+{\mathrm{B}}^{+--}+{\mathrm{C}}^{-+-}+{\mathrm{D}}^{--+}\\ +{\mathrm{E}}^{++-}+{\mathrm{F}}^{+-+}+{\mathrm{G}}^{-++}+{\mathrm{H}}^{+++}\end{array}\right|}^2+\\ 4{\mu }_{\mathrm{c}}{\left|{\mathrm{k}}_{\mathrm{x}}{\mathrm{k}}_{\mathrm{y}}\right|}^2{\left|\begin{array}{c}{\mathrm{A}}^{---}-{\mathrm{B}}^{+--}-{\mathrm{C}}^{-+-}+{\mathrm{D}}^{--+}\\ +{\mathrm{E}}^{++-}-{\mathrm{F}}^{+-+}-{\mathrm{G}}^{-++}+{\mathrm{H}}^{+++}\end{array}\right|}^2+\\ 4{\mu }_{\mathrm{c}}{\left|{\mathrm{k}}_{\mathrm{y}}{\mathrm{k}}_{\mathrm{z}}\right|}^2{\left|\begin{array}{c}{\mathrm{A}}^{---}+{\mathrm{B}}^{+--}-{\mathrm{C}}^{-+-}-{\mathrm{D}}^{--+}\\ -{\mathrm{E}}^{++-}-{\mathrm{F}}^{+-+}+{\mathrm{G}}^{-++}+{\mathrm{H}}^{+++}\end{array}\right|}^2+\\ 4{\mu }_{\mathrm{c}}{\left|{\mathrm{k}}_{\mathrm{z}}{\mathrm{k}}_{\mathrm{x}}\right|}^2{\left|\begin{array}{c}{\mathrm{A}}^{---}-{\mathrm{B}}^{+--}+{\mathrm{C}}^{-+-}-{\mathrm{D}}^{--+}\\ -{\mathrm{E}}^{++-}+{\mathrm{F}}^{+-+}-{\mathrm{G}}^{-++}+{\mathrm{H}}^{+++}\end{array}\right|}^2\end{array}\right]$$

(17)

$$〈{\mathrm{e}}_{\mathrm{k}}〉=\frac{1}{4}{\mathrm{\rho \omega }}^2\mathrm{Re}\left[\begin{array}{l}{\left|{\mathrm{k}}_{\mathrm{x}}\right|}^2{\left|\begin{array}{c}{\mathrm{A}}^{---}-{\mathrm{B}}^{+--}+{\mathrm{C}}^{-+-}+{\mathrm{D}}^{--+}\\ -{\mathrm{E}}^{++-}-{\mathrm{F}}^{+-+}+{\mathrm{G}}^{-++}-{\mathrm{H}}^{+++}\end{array}\right|}^2+\\ {\left|{\mathrm{k}}_{\mathrm{y}}\right|}^2{\left|\begin{array}{c}{\mathrm{A}}^{---}+{\mathrm{B}}^{+--}-{\mathrm{C}}^{-+-}+{\mathrm{D}}^{--+}\\ -{\mathrm{E}}^{++-}+{\mathrm{F}}^{+-+}-{\mathrm{G}}^{-++}-{\mathrm{H}}^{+++}\end{array}\right|}^2+\\ {\left|{\mathrm{k}}_{\mathrm{z}}\right|}^2{\left|\begin{array}{c}{\mathrm{A}}^{---}+{\mathrm{B}}^{+--}+{\mathrm{C}}^{-+-}-{\mathrm{D}}^{--+}\\ +{\mathrm{E}}^{++-}-{\mathrm{F}}^{+-+}-{\mathrm{G}}^{-++}-{\mathrm{H}}^{+++}\end{array}\right|}^2\end{array}\right]$$

(18)

For the elastic solid medium, the total energy density is the linear combination of the kinetic and potential energy densities. The time-averaged total energy density of the solid can be expressed as

$$〈{\mathrm{e}}_{\mathrm{t}}〉=\frac{1}{4}\mathrm{Re}\left[\begin{array}{l}\left\{{\mathrm{\lambda }}_{\mathrm{c}}{\left|{\mathrm{k}}_{\mathrm{L}}^2\right|}^2+2{\mu }_{\mathrm{c}}\left({\left|{\mathrm{k}}_{\mathrm{x}}^2\right|}^2+{\left|{\mathrm{k}}_{\mathrm{y}}^2\right|}^2+{\left|{\mathrm{k}}_{\mathrm{z}}^2\right|}^2\right)\right\}\\ \times {\left|\begin{array}{c}{\mathrm{A}}^{---}+{\mathrm{B}}^{+--}+{\mathrm{C}}^{-+-}+{\mathrm{D}}^{--+}\\ +{\mathrm{E}}^{++-}+{\mathrm{F}}^{+-+}+{\mathrm{G}}^{-++}+{\mathrm{H}}^{+++}\end{array}\right|}^2\\ +4{\mu }_{\mathrm{c}}{\left|{\mathrm{k}}_{\mathrm{x}}{\mathrm{k}}_{\mathrm{y}}\right|}^2{\left|\begin{array}{c}{\mathrm{A}}^{---}-{\mathrm{B}}^{+--}-{\mathrm{C}}^{-+-}+{\mathrm{D}}^{--+}\\ +{\mathrm{E}}^{++-}-{\mathrm{F}}^{+-+}-{\mathrm{G}}^{-++}+{\mathrm{H}}^{+++}\end{array}\right|}^2\\ +4{\mu }_{\mathrm{c}}{\left|{\mathrm{k}}_{\mathrm{y}}{\mathrm{k}}_{\mathrm{z}}\right|}^2{\left|\begin{array}{c}{\mathrm{A}}^{---}+{\mathrm{B}}^{+--}-{\mathrm{C}}^{-+-}-{\mathrm{D}}^{--+}\\ -{\mathrm{E}}^{++-}-{\mathrm{F}}^{+-+}+{\mathrm{G}}^{-++}+{\mathrm{H}}^{+++}\end{array}\right|}^2\\ +4{\mu }_{\mathrm{c}}{\left|{\mathrm{k}}_{\mathrm{z}}{\mathrm{k}}_{\mathrm{x}}\right|}^2{\left|\begin{array}{c}{\mathrm{A}}^{---}-{\mathrm{B}}^{+--}+{\mathrm{C}}^{-+-}-{\mathrm{D}}^{--+}\\ -{\mathrm{E}}^{++-}+{\mathrm{F}}^{+-+}-{\mathrm{G}}^{-++}+{\mathrm{H}}^{+++}\end{array}\right|}^2\\ +{\mathrm{\rho \omega }}^2{\left|{\mathrm{k}}_{\mathrm{x}}\right|}^2{\left|\begin{array}{c}{\mathrm{A}}^{---}-{\mathrm{B}}^{+--}+{\mathrm{C}}^{-+-}+{\mathrm{D}}^{--+}\\ -{\mathrm{E}}^{++-}-{\mathrm{F}}^{+-+}+{\mathrm{G}}^{-++}-{\mathrm{H}}^{+++}\end{array}\right|}^2\\ +{\mathrm{\rho \omega }}^2{\left|{\mathrm{k}}_{\mathrm{y}}\right|}^2{\left|\begin{array}{c}{\mathrm{A}}^{---}+{\mathrm{B}}^{+--}-{\mathrm{C}}^{-+-}+{\mathrm{D}}^{--+}\\ -{\mathrm{E}}^{++-}+{\mathrm{F}}^{+-+}-{\mathrm{G}}^{-++}-{\mathrm{H}}^{+++}\end{array}\right|}^2\\ +{\mathrm{\rho \omega }}^2{\left|{\mathrm{k}}_{\mathrm{z}}\right|}^2{\left|\begin{array}{c}{\mathrm{A}}^{---}+{\mathrm{B}}^{+--}+{\mathrm{C}}^{-+-}-{\mathrm{D}}^{--+}\\ +{\mathrm{E}}^{++-}-{\mathrm{F}}^{+-+}-{\mathrm{G}}^{-++}-{\mathrm{H}}^{+++}\end{array}\right|}^2\end{array}\right]$$

(19)

Substituting equation (13) into equations (10)–(12), the expanded expressions of each components of the time-averaged intensity are derived as

$$〈{\mathrm{I}}_{\mathrm{x}}〉=\frac{\mathrm{\omega }}{2}\mathrm{Re}\left[\begin{array}{l}\left({\mathrm{\lambda }}_{\mathrm{c}}{\mathrm{k}}_{\mathrm{L}}^2{{\mathrm{k}}_{\mathrm{x}}}^{\mathrm{*}}+2{\mu }_{\mathrm{c}}{\mathrm{k}}_{\mathrm{x}}{\left|{\mathrm{k}}_{\mathrm{x}}\right|}^2\right)\\ \times \left(\begin{array}{c}{\mathrm{A}}^{---}+{\mathrm{B}}^{+--}+{\mathrm{C}}^{-+-}+{\mathrm{D}}^{--+}\\ +{\mathrm{E}}^{++-}+{\mathrm{F}}^{+-+}+{\mathrm{G}}^{-++}+{\mathrm{H}}^{+++}\end{array}\right)\\ \times {\left(\begin{array}{c}{\mathrm{A}}^{---}-{\mathrm{B}}^{+--}+{\mathrm{C}}^{-+-}+{\mathrm{D}}^{--+}\\ -{\mathrm{E}}^{++-}-{\mathrm{F}}^{+-+}+{\mathrm{G}}^{-++}-{\mathrm{H}}^{+++}\end{array}\right)}^{\mathrm{*}}+\\ 2{\mu }_{\mathrm{c}}{\mathrm{k}}_{\mathrm{x}}{\left|{\mathrm{k}}_{\mathrm{y}}\right|}^2\left(\begin{array}{c}{\mathrm{A}}^{---}-{\mathrm{B}}^{+--}-{\mathrm{C}}^{-+-}+{\mathrm{D}}^{--+}\\ +{\mathrm{E}}^{++-}-{\mathrm{F}}^{+-+}-{\mathrm{G}}^{-++}+{\mathrm{H}}^{+++}\end{array}\right)\\ \times {\left(\begin{array}{c}{\mathrm{A}}^{---}+{\mathrm{B}}^{+--}-{\mathrm{C}}^{-+-}+{\mathrm{D}}^{--+}\\ -{\mathrm{E}}^{++-}+{\mathrm{F}}^{+-+}-{\mathrm{G}}^{-++}-{\mathrm{H}}^{+++}\end{array}\right)}^{\mathrm{*}}+\\ 2{\mu }_{\mathrm{c}}{\mathrm{k}}_{\mathrm{x}}{\left|{\mathrm{k}}_{\mathrm{z}}\right|}^2\left(\begin{array}{c}{\mathrm{A}}^{---}-{\mathrm{B}}^{+--}+{\mathrm{C}}^{-+-}-{\mathrm{D}}^{--+}\\ -{\mathrm{E}}^{++-}+{\mathrm{F}}^{+-+}-{\mathrm{G}}^{-++}+{\mathrm{H}}^{+++}\end{array}\right)\\ \times {\left(\begin{array}{c}{\mathrm{A}}^{---}+{\mathrm{B}}^{+--}+{\mathrm{C}}^{-+-}-{\mathrm{D}}^{--+}\\ +{\mathrm{E}}^{++-}-{\mathrm{F}}^{+-+}-{\mathrm{G}}^{-++}-{\mathrm{H}}^{+++}\end{array}\right)}^{\mathrm{*}}\end{array}\right]$$

(20)

$$〈{\mathrm{I}}_{\mathrm{y}}〉=\frac{\mathrm{\omega }}{2}\mathrm{Re}\left[\begin{array}{l}\left({\mathrm{\lambda }}_{\mathrm{c}}{\mathrm{k}}_{\mathrm{L}}^2{{\mathrm{k}}_{\mathrm{y}}}^{\mathrm{*}}+2{\mu }_{\mathrm{c}}{\mathrm{k}}_{\mathrm{y}}{\left|{\mathrm{k}}_{\mathrm{y}}\right|}^2\right)\\ \times \left(\begin{array}{c}{\mathrm{A}}^{---}+{\mathrm{B}}^{+--}+{\mathrm{C}}^{-+-}+{\mathrm{D}}^{--+}\\ +{\mathrm{E}}^{++-}+{\mathrm{F}}^{+-+}+{\mathrm{G}}^{-++}+{\mathrm{H}}^{+++}\end{array}\right)\\ \times {\left(\begin{array}{c}{\mathrm{A}}^{---}+{\mathrm{B}}^{+--}-{\mathrm{C}}^{-+-}+{\mathrm{D}}^{--+}\\ -{\mathrm{E}}^{++-}+{\mathrm{F}}^{+-+}-{\mathrm{G}}^{-++}-{\mathrm{H}}^{+++}\end{array}\right)}^{\mathrm{*}}+\\ 2{\mu }_{\mathrm{c}}{\mathrm{k}}_{\mathrm{y}}{\left|{\mathrm{k}}_{\mathrm{x}}\right|}^2\left(\begin{array}{c}{\mathrm{A}}^{---}-{\mathrm{B}}^{+--}-{\mathrm{C}}^{-+-}+{\mathrm{D}}^{--+}\\ +{\mathrm{E}}^{++-}-{\mathrm{F}}^{+-+}-{\mathrm{G}}^{-++}+{\mathrm{H}}^{+++}\end{array}\right)\\ \times {\left(\begin{array}{c}{\mathrm{A}}^{---}+{\mathrm{B}}^{+--}-{\mathrm{C}}^{-+-}+{\mathrm{D}}^{--+}\\ -{\mathrm{E}}^{++-}+{\mathrm{F}}^{+-+}-{\mathrm{G}}^{-++}-{\mathrm{H}}^{+++}\end{array}\right)}^{\mathrm{*}}+\\ 2{\mu }_{\mathrm{c}}{\mathrm{k}}_{\mathrm{y}}{\left|{\mathrm{k}}_{\mathrm{z}}\right|}^2\left(\begin{array}{c}{\mathrm{A}}^{---}+{\mathrm{B}}^{+--}-{\mathrm{C}}^{-+-}-{\mathrm{D}}^{--+}\\ -{\mathrm{E}}^{++-}-{\mathrm{F}}^{+-+}+{\mathrm{G}}^{-++}+{\mathrm{H}}^{+++}\end{array}\right)\\ \times {\left(\begin{array}{c}{\mathrm{A}}^{---}+{\mathrm{B}}^{+--}+{\mathrm{C}}^{-+-}-{\mathrm{D}}^{--+}\\ +{\mathrm{E}}^{++-}-{\mathrm{F}}^{+-+}-{\mathrm{G}}^{-++}-{\mathrm{H}}^{+++}\end{array}\right)}^{\mathrm{*}}\end{array}\right]$$

(21)

$$〈{\mathrm{I}}_{\mathrm{z}}〉=\frac{\mathrm{\omega }}{2}\mathrm{Re}\left[\begin{array}{l}\left({\mathrm{\lambda }}_{\mathrm{c}}{\mathrm{k}}_{\mathrm{L}}^2{{\mathrm{k}}_{\mathrm{z}}}^{\mathrm{*}}+2{\mu }_{\mathrm{c}}{\mathrm{k}}_{\mathrm{z}}{\left|{\mathrm{k}}_{\mathrm{z}}\right|}^2\right)\\ \times \left(\begin{array}{c}{\mathrm{A}}^{---}+{\mathrm{B}}^{+--}+{\mathrm{C}}^{-+-}+{\mathrm{D}}^{--+}\\ +{\mathrm{E}}^{++-}+{\mathrm{F}}^{+-+}+{\mathrm{G}}^{-++}+{\mathrm{H}}^{+++}\end{array}\right)\\ \times {\left(\begin{array}{c}{\mathrm{A}}^{---}+{\mathrm{B}}^{+--}+{\mathrm{C}}^{-+-}-{\mathrm{D}}^{--+}\\ +{\mathrm{E}}^{++-}-{\mathrm{F}}^{+-+}-{\mathrm{G}}^{-++}-{\mathrm{H}}^{+++}\end{array}\right)}^{\mathrm{*}}+\\ 2{\mu }_{\mathrm{c}}{\mathrm{k}}_{\mathrm{z}}{\left|{\mathrm{k}}_{\mathrm{x}}\right|}^2\left(\begin{array}{c}{\mathrm{A}}^{---}-{\mathrm{B}}^{+--}+{\mathrm{C}}^{-+-}-{\mathrm{D}}^{--+}\\ -{\mathrm{E}}^{++-}+{\mathrm{F}}^{+-+}-{\mathrm{G}}^{-++}+{\mathrm{H}}^{+++}\end{array}\right)\\ \times {\left(\begin{array}{c}{\mathrm{A}}^{---}-{\mathrm{B}}^{+--}+{\mathrm{C}}^{-+-}+{\mathrm{D}}^{--+}\\ -{\mathrm{E}}^{++-}-{\mathrm{F}}^{+-+}+{\mathrm{G}}^{-++}-{\mathrm{H}}^{+++}\end{array}\right)}^{\mathrm{*}}+\\ 2{\mu }_{\mathrm{c}}{\mathrm{k}}_{\mathrm{z}}{\left|{\mathrm{k}}_{\mathrm{y}}\right|}^2\left(\begin{array}{c}{\mathrm{A}}^{---}+{\mathrm{B}}^{+--}-{\mathrm{C}}^{-+-}-{\mathrm{D}}^{--+}\\ -{\mathrm{E}}^{++-}-{\mathrm{F}}^{+-+}+{\mathrm{G}}^{-++}+{\mathrm{H}}^{+++}\end{array}\right)\\ \times {\left(\begin{array}{c}{\mathrm{A}}^{---}+{\mathrm{B}}^{+--}-{\mathrm{C}}^{-+-}+{\mathrm{D}}^{--+}\\ -{\mathrm{E}}^{++-}+{\mathrm{F}}^{+-+}-{\mathrm{G}}^{-++}-{\mathrm{H}}^{+++}\end{array}\right)}^{\mathrm{*}}\end{array}\right]$$

(22)

The relation equation between the energy density and intensity is required to obtain the energy governing equation. Here, no obvious relations can be found between the energy density and the intensity. The energy and intensity expressions can be simplified using a two-step process. The first is to neglect all the terms of order ${\eta }^2$ and higher, which is proper if the hysteretic damping for the solid is sufficiently small. The second is to neglect all terms containing sinusoidal functions of the wavenumber. This simplification is the equivalent to applying space-averaging of the form

$$〈\overline{{\mathrm{e}}_{\mathrm{t}}}〉=\frac{{\mathrm{k}}_{\mathrm{x}1}{\mathrm{k}}_{\mathrm{y}1}{\mathrm{k}}_{\mathrm{z}1}}{{\pi }^3}\underset{0}{\overset{\pi /{\mathrm{k}}_{\mathrm{z}1}}{\int }}\underset{0}{\overset{\pi /{\mathrm{k}}_{\mathrm{y}1}}{\int }}\underset{0}{\overset{\pi /{\mathrm{k}}_{\mathrm{x}1}}{\int }}〈{\mathrm{e}}_{\mathrm{t}}〉\mathrm{dx}\mathrm{dy}\mathrm{dz}$$

(23)

$$〈\overline{{\mathrm{I}}_{\mathrm{i}}}〉=\frac{{\mathrm{k}}_{\mathrm{x}1}{\mathrm{k}}_{\mathrm{y}1}{\mathrm{k}}_{\mathrm{z}1}}{{\pi }^3}\underset{0}{\overset{\pi /{\mathrm{k}}_{\mathrm{z}1}}{\int }}\underset{0}{\overset{\pi /{\mathrm{k}}_{\mathrm{y}1}}{\int }}\underset{0}{\overset{\pi /{\mathrm{k}}_{\mathrm{x}1}}{\int }}〈{\mathrm{I}}_{\mathrm{i}}〉\mathrm{dx}\mathrm{dy}\mathrm{dz}$$

(24)

where $\overline{〈{\mathrm{e}}_{\mathrm{t}}〉}$ and $\overline{〈{\mathrm{I}}_{\mathrm{i}}〉}$ are the time- and space-averaged energy density and Cartesian component of intensity, and subscript $\mathrm{i}$ indicates the x, y, z component of the intensity, respectively. Simplified expressions for the time- and space-averaged energy density and intensity are written as

$$〈\overline{{\mathrm{e}}_{\mathrm{t}}}〉=\frac{1}{2}(\mathrm{\lambda }+2\mu ){\mathrm{k}}_{\mathrm{L}1}^4\times \left(\begin{array}{c}{\left|\mathrm{A}\right|}^2{\mathrm{e}}^{---}+{\left|\mathrm{B}\right|}^2{\mathrm{e}}^{+--}+{\left|\mathrm{C}\right|}^2{\mathrm{e}}^{-+-}+{\left|\mathrm{D}\right|}^2{\mathrm{e}}^{--+}\\ +{\left|\mathrm{E}\right|}^2{\mathrm{e}}^{++-}+{\left|\mathrm{F}\right|}^2{\mathrm{e}}^{+-+}+{\left|\mathrm{G}\right|}^2{\mathrm{e}}^{-++}+{\left|\mathrm{H}\right|}^2{\mathrm{e}}^{+++}\end{array}\right)$$

(25)

$$〈\overline{{\mathrm{I}}_{\mathrm{x}}}〉=\frac{1}{2}{\mathrm{\omega k}}_{\mathrm{x}1}(\mathrm{\lambda }+2\mu ){\mathrm{k}}_{\mathrm{L}1}^2\times \left(\begin{array}{c}{\left|\mathrm{A}\right|}^2{\mathrm{e}}^{---}-{\left|\mathrm{B}\right|}^2{\mathrm{e}}^{+--}+{\left|\mathrm{C}\right|}^2{\mathrm{e}}^{-+-}+{\left|\mathrm{D}\right|}^2{\mathrm{e}}^{--+}\\ -{\left|\mathrm{E}\right|}^2{\mathrm{e}}^{++-}-{\left|\mathrm{F}\right|}^2{\mathrm{e}}^{+-+}+{\left|\mathrm{G}\right|}^2{\mathrm{e}}^{-++}-{\left|\mathrm{H}\right|}^2{\mathrm{e}}^{+++}\end{array}\right)$$

(26)

$$〈\overline{{\mathrm{I}}_{\mathrm{y}}}〉=\frac{1}{2}{\mathrm{\omega k}}_{\mathrm{y}1}(\mathrm{\lambda }+2\mu ){\mathrm{k}}_{\mathrm{L}1}^2\times \left(\begin{array}{c}{\left|\mathrm{A}\right|}^2{\mathrm{e}}^{---}+{\left|\mathrm{B}\right|}^2{\mathrm{e}}^{+--}-{\left|\mathrm{C}\right|}^2{\mathrm{e}}^{-+-}+{\left|\mathrm{D}\right|}^2{\mathrm{e}}^{--+}\\ -{\left|\mathrm{E}\right|}^2{\mathrm{e}}^{++-}+{\left|\mathrm{F}\right|}^2{\mathrm{e}}^{+-+}-{\left|\mathrm{G}\right|}^2{\mathrm{e}}^{-++}-{\left|\mathrm{H}\right|}^2{\mathrm{e}}^{+++}\end{array}\right)$$

(27)

$$〈\overline{{\mathrm{I}}_{\mathrm{z}}}〉=\frac{1}{2}{\mathrm{\omega k}}_{\mathrm{z}1}(\mathrm{\lambda }+2\mu ){\mathrm{k}}_{\mathrm{L}1}^2\times \left(\begin{array}{c}{\left|\mathrm{A}\right|}^2{\mathrm{e}}^{---}+{\left|\mathrm{B}\right|}^2{\mathrm{e}}^{+--}+{\left|\mathrm{C}\right|}^2{\mathrm{e}}^{-+-}-{\left|\mathrm{D}\right|}^2{\mathrm{e}}^{--+}\\ +{\left|\mathrm{E}\right|}^2{\mathrm{e}}^{++-}-{\left|\mathrm{F}\right|}^2{\mathrm{e}}^{+-+}-{\left|\mathrm{G}\right|}^2{\mathrm{e}}^{-++}-{\left|\mathrm{H}\right|}^2{\mathrm{e}}^{+++}\end{array}\right)$$

(28)

where $e^{\pm \pm \pm }$ are defined as $\exp (\pm \eta {\mathrm{k}}_{\mathrm{x}1}\mathrm{x}\pm \eta {\mathrm{k}}_{\mathrm{y}1}\mathrm{y}\pm \eta {\mathrm{k}}_{\mathrm{z}1}\mathrm{z})$ . These simplified expressions are relatively simple functions of A, B, …, H, and this reveals that the time- and space-averaged energy density and intensity components are related in the following manner

$$〈\overline{{\mathrm{I}}_{\mathrm{x}}}〉=-\frac{{\mathrm{c}}_{\mathrm{L}}^2}{\eta \omega }\frac{\partial 〈\overline{{\mathrm{e}}_{\mathrm{t}}}〉}{\partial \mathrm{x}}$$

(29)

$$〈\overline{{\mathrm{I}}_{\mathrm{y}}}〉=-\frac{{\mathrm{c}}_{\mathrm{L}}^2}{\eta \omega }\frac{\partial 〈\overline{{\mathrm{e}}_{\mathrm{t}}}〉}{\partial \mathrm{y}}$$

(30)

$$〈\overline{{\mathrm{I}}_{\mathrm{z}}}〉=-\frac{{\mathrm{c}}_{\mathrm{L}}^2}{\eta \omega }\frac{\partial 〈\overline{{\mathrm{e}}_{\mathrm{t}}}〉}{\partial \mathrm{z}}$$

(31)

Consequently, rearranging equations (29)–(31), the relation equation between the time- and space-averaged energy density and intensity yields

$$\overline{〈\overrightarrow{\mathrm{I}}〉}=-\frac{{\mathrm{c}}_{\mathrm{L}}^2}{\eta \omega }\nabla 〈\overline{{\mathrm{e}}_{\mathrm{t}}}〉$$

(32)

For a steady-state condition, the energy balance equation of an elastic solid medium is expressed as

$$\nabla ·\overrightarrow{\mathrm{I}}+{\pi }_{\mathrm{diss}}={\pi }_{\mathrm{in}}$$

(33)

where ${\pi }_{\mathrm{diss}}$ and ${\pi }_{\mathrm{in}}$ are the time-averaged dissipated power at a point due to structural damping of the system and time-averaged input power, respectively. In the case of light hysteretic damping, the dissipated power at medium vibrating with a circular frequency, $\omega$ , is proportional to the structural damping loss factor, and the reversible energy density in the medium in such a way that

$${\pi }_{\mathrm{diss}}=\eta \omega 〈\mathrm{e}〉$$

(34)

The gradient of the intensity in equation (33) can be expressed in terms of the second derivative of the energy density using the relation equation between the intensity and the energy density in equation (32), such that the energy governing equation for the dilatational wave in the elastic solid medium without input power can be derived as

$$-\frac{{\mathrm{c}}_{\mathrm{gL}}^2}{\eta \omega }{\nabla }^2〈\overline{{\mathrm{e}}_{\mathrm{t}}}〉+\eta \omega 〈\overline{{\mathrm{e}}_{\mathrm{t}}}〉=0$$

(35)

Energy flow solutions for dilatational waves in the elastic solid medium

This section performs an analysis of a vibrating solid cube excited by spherical cavity source subjected to internal pressure by solving equation (35) using a Fourier series technique to find the energy density distributions. Since this study develops energy flow model for the dilatational wave in 3D elastic solids, the conditions of the structure are set so that only the dilatational wave exists in the structure, not the distortion wave.

The model used in these analyses is cube-shaped and contains a spherical cavity source at position $\mathrm{x}={\mathrm{x}}_0$ , $\mathrm{y}={\mathrm{y}}_0$ , $\mathrm{z}={\mathrm{z}}_0$ , constrained by a traction-free boundary condition on all six faces of the cube (Figure 1). For the traction-free boundary condition that corresponds to an elastic half-space constrained by a rigid, lubricated boundary, mode conversion of the wave type does not occur. ¹⁶ This lack of occurrence of mode conversion implies that the wave type of the reflected wave equals that of the incident wave.

Figure 1.

Cube shape structure with spherical cavity.

The spherical symmetric cavity surface excited by time-varying-but-uniformly applied normal pressure emanates only the dilatational wave, and the six faces of the cube are constrained by the traction-free condition, so the resulting waves generated from the spherical source and reflecting waves from the boundary can exist only as dilatational waves. Thus, it is reasonable to assume that there is no distortion wave for this example. If the diameter of the source is relatively smaller than the cube dimension, the spherical cavity source can be considered as simple monopole source.

When the spherical source is applied at $\mathrm{x}={\mathrm{x}}_0$ , $\mathrm{y}={\mathrm{y}}_0$ , $\mathrm{z}={\mathrm{z}}_0$ of the solid, the energy governing equation of the vibrating cube becomes

$$-\frac{{\mathrm{c}}_{\mathrm{gL}}^2}{\eta \omega }{\nabla }^2〈\overline{{\mathrm{e}}_{\mathrm{t}}}〉+\eta \omega 〈\overline{{\mathrm{e}}_{\mathrm{t}}}〉={\Pi }_{\mathrm{in}}\delta (\mathrm{x}-{\mathrm{x}}_0)\delta (\mathrm{y}-{\mathrm{y}}_0)\delta (\mathrm{z}-{\mathrm{z}}_0)$$

(36)

where ${\Pi }_{\mathrm{in}}$ is the time-averaged input power due to the internal pressure, and $\delta$ is the Dirac delta. On boundary surface, which is subjected to the traction-free boundary condition, it can be assumed that there is no normal component for the intensity, and such condition is equivalent to that in which the face normal component of the partial derivative of the energy density is zero. To satisfy the boundary condition, the general solution of equation (35), the energy density, is expressed using the triple Fourier series of the cosine functions with respect to the spatial variables, $\mathrm{x}$ , $\mathrm{y}$ , $\mathrm{z}$

$$〈\overline{{\mathrm{e}}_{\mathrm{t}}}〉=\sum _{\mathrm{l}=0}^{\infty }\sum _{\mathrm{m}=0}^{\infty }\sum _{\mathrm{n}=0}^{\infty }{\mathrm{g}}_{\mathrm{nml}}\cos \left(\frac{\mathrm{n}\mathrm{\pi }\mathrm{x}}{{\mathrm{L}}_{\mathrm{x}}}\right)\cos \left(\frac{\mathrm{m}\mathrm{\pi }\mathrm{y}}{{\mathrm{L}}_{\mathrm{y}}}\right)\cos \left(\frac{\mathrm{l}\mathrm{\pi }\mathrm{z}}{{\mathrm{L}}_{\mathrm{z}}}\right)$$

(37)

where ${\mathrm{g}}_{\mathrm{nml}}$ is the $(\mathrm{n},\mathrm{m},\mathrm{l})$ mode coefficient of the energy density solution, and ${\mathrm{L}}_{\mathrm{x}}$ , ${\mathrm{L}}_{\mathrm{y}}$ , and ${\mathrm{L}}_{\mathrm{z}}$ are the dimensions of the cube. The time-averaged input power, ${\Pi }_{\mathrm{in}}$ is estimated as

$${\Pi }_{\mathrm{in}}=\frac{1}{2}\mathrm{Re}\left[4\pi {{\mathrm{r}}_0}^2\mathrm{p}{\frac{\partial {{\mathrm{u}}_{\mathrm{r}}}^{*}}{\partial \mathrm{t}}}_{\mathrm{at}({\mathrm{x}}_0,{\mathrm{y}}_0,{\mathrm{z}}_0)}\right]$$

(38)

where ${\mathrm{r}}_0$ , $\mathrm{p}$ , and ${\mathrm{u}}_{\mathrm{r}}$ are the radius of the spherical cavity source, harmonic internal pressure on the spherical cavity, and radial oscillatory surface displacement at the cavity surface, respectively. Substituting equations (38) and (37) into equation (36) yields

$${\mathrm{g}}_{\mathrm{nml}}=\frac{\eta \omega }{{\mathrm{c}}_{\mathrm{gL}}^2}\frac{{\mathrm{G}}_{\mathrm{nml}}}{[{(\mathrm{n}\mathrm{\pi }/{\mathrm{L}}_{\mathrm{x}})}^2+{(\mathrm{m}\mathrm{\pi }/{\mathrm{L}}_{\mathrm{y}})}^2+{(\mathrm{l}\mathrm{\pi }/{\mathrm{L}}_{\mathrm{z}})}^2+{(\eta \omega /{\mathrm{c}}_{\mathrm{gL}})}^2]}$$

(39)

where ${\mathrm{G}}_{\mathrm{nml}}$ and $\mathrm{k}$ are

$${\mathrm{G}}_{\mathrm{nml}}=\frac{\mathrm{k}{\Pi }_{\mathrm{in}}}{{\mathrm{L}}_{\mathrm{x}}{\mathrm{L}}_{\mathrm{y}}{\mathrm{L}}_{\mathrm{z}}}\cos \left(\frac{\mathrm{n}\mathrm{\pi }{\mathrm{x}}_0}{{\mathrm{L}}_{\mathrm{x}}}\right)\cos \left(\frac{\mathrm{m}\mathrm{\pi }{\mathrm{y}}_0}{{\mathrm{L}}_{\mathrm{y}}}\right)\cos \left(\frac{\mathrm{l}\mathrm{\pi }{\mathrm{z}}_0}{{\mathrm{L}}_{\mathrm{z}}}\right)$$

(40)

$$\mathrm{k}={\epsilon }_{\mathrm{n}}{\epsilon }_{\mathrm{m}}{\epsilon }_{\mathrm{l}},{\epsilon }_{\mathrm{i}}=\{\begin{array}{c}1(i=0)\hfill \\ \multicolumn{1}{c}{2(i\ne 0)}\end{array}$$

(41)

Validation

In this section, the results of the energy density obtained from the EFA are compared to those calculated using NASTRAN. In general, the FEA is optimized at the low-frequency range because the lower the frequency, the smaller the error. Contrary, because of the statistical approach, the EFA has the characteristic that the error decreases as the frequency increases. Because of the different characteristics of the two methods, in order to validate EFA with FEA, verification frequency was selected as the highest frequency of the reliable frequency range of FEA.

The cube model that is applied in this verification has a size of $200\mathrm{m}\times 200\mathrm{m}\times 200\mathrm{m}$ , and it contains a spherical cavity in its center. The radius of the spherical cavity is ${\mathrm{r}}_0=2\mathrm{m}$ , and the location of the spherical cavity source at the center can be written as

$${\mathrm{x}}_0=\frac{{\mathrm{L}}_{\mathrm{x}}}{2},{\mathrm{y}}_0=\frac{{\mathrm{L}}_{\mathrm{y}}}{2},{\mathrm{z}}_0=\frac{{\mathrm{L}}_{\mathrm{z}}}{2}$$

(42)

The material properties are set to those of steel (Young’s modulus $\mathrm{E}=2.07\times {10}^{11}\mathrm{N}/{\mathrm{m}}^2$ , density $\rho =7800\mathrm{kg}/{\mathrm{m}}^3$ , and Poisson’s ratio $\nu =0.28$ ), and the amplitude of the harmonic internal pressure $\mathrm{p}$ on the spherical cavity is ${10}^4\mathrm{N}/{\mathrm{m}}^2$ .

Figure 2 shows the FE model used for NASTRAN analysis. Since the shape and boundary conditions of the structure are symmetrical, one-eighth of the cube is modeled to improve the efficiency of the analysis, and the convergence test was performed according to the number of elements. The total number of elements constituting the FE model used for the verification is about 420,000.

Figure 2.

FE model used for NASTRAN analysis.

In the first example, the exciting frequency is $\mathrm{f}=500\mathrm{Hz}$ , and the structural damping coefficient is $\eta =0.06$ . The spatial distribution of the energy density obtained from the energy flow model and NASTRAN are shown in Figures 3 and 4. Figure 3 shows the energy density distribution on a quarter of the mid-cross section, and Figure 4 shows the distribution of the energy density where the observation points stand in the line that connects the center point of the cube to the center point of a surface of the cube. The energy density of the energy flow model, as shown in Figure 3(a), varied smoothly in space, without any fluctuation. The energy density obtained from NASTRAN decreased universally as the distance from the excitation location increased, and it fluctuated locally in space, as shown in Figure 3(b). A comparison of the results of the energy flow model and NASTRAN reveals that the energy flow model solution represents the global variation of the NASTRAN solution quite well. A comparison with the qualified results from NASTRAN validates the use of the energy flow model for the 3D elastic solid.

Figure 3.

Energy density distribution of the elastic cube on quarter of mid-cross section when $\mathrm{f}=500\mathrm{Hz}$ and $\mathrm{\eta }=0.06$ . The reference energy density is $1\times {10}^{-12}\mathrm{J}/{\mathrm{m}}^3$ . ${\mathrm{L}}_{\mathrm{x}}={\mathrm{L}}_{\mathrm{y}}={\mathrm{L}}_{\mathrm{z}}=200\mathrm{m}$ : (a) energy flow solutions and (b) NASTRAN solution.

Figure 4.

Energy density distribution of the elastic cube when $\mathrm{f}=500\mathrm{Hz}$ and $\mathrm{\eta }=0.06$ . The reference energy density is $1\times {10}^{-12}\mathrm{J}/{\mathrm{m}}^3$ . ${\mathrm{L}}_{\mathrm{x}}={\mathrm{L}}_{\mathrm{y}}={\mathrm{L}}_{\mathrm{z}}=200\mathrm{m}$ ; ——: EFA solution; ○: NASTRAN solution.

FEA is inapplicable in high-frequency problem because the uncertainty of the system increases and the computational error is accumulated as the frequency increases. However, since the EFA is based on statistical assumptions, the error is reduced as the frequency increases. Therefore, it is possible to analyze the medium- to high-frequency range that cannot be performed by conventional FEA or BEA.

The results obtained from the energy flow method for the high-frequency problem are shown below. In the second example, EFAs at two different frequencies $\mathrm{f}=700\mathrm{Hz}$ and $1000\mathrm{Hz}$ are implemented.

The distributions of the energy densities on a quarter of the mid-section are shown in Figures 5 and 6. The results of the energy densities calculated from the EFA do not fluctuate locally, but decrease smoothly throughout the whole structure, as expected. The maximum/minimum values of the energy density at frequencies of $\mathrm{f}=700\mathrm{Hz}$ and $1000\mathrm{Hz}$ are $64.0/44.1\mathrm{dB}$ and $66.6/39.3\mathrm{dB}$ , respectively. A comparison of the two results of this example shows that the energy difference between the maximum and minimum values increases as the frequency increases, which implies that the energy density spatially decreases more quickly at a higher frequency.

Figure 5.

Energy density distribution of the elastic cube on quarter of mid-cross section when $\mathrm{f}=700\mathrm{Hz}$ and $\mathrm{\eta }=0.03$ . The reference energy density is $1\times {10}^{-12}\mathrm{J}/{\mathrm{m}}^3$ . ${\mathrm{L}}_{\mathrm{x}}={\mathrm{L}}_{\mathrm{y}}={\mathrm{L}}_{\mathrm{z}}=200\mathrm{m}$ .

Figure 6.

Energy density distribution of the elastic cube on quarter of mid-cross section when $\mathrm{f}=1000\mathrm{Hz}$ and $\eta =0.03$ . The reference energy density is $1\times {10}^{-12}\mathrm{J}/{\mathrm{m}}^3$ . ${\mathrm{L}}_{\mathrm{x}}={\mathrm{L}}_{\mathrm{y}}={\mathrm{L}}_{\mathrm{z}}=200\mathrm{m}$ .

The last example implements EFAs for two different damping factors. Figures 7 and 8 show the results of the energy density for damping loss factors eta = 0.01 and 0.06, respectively, at a frequency f = 700 Hz. The maximum/minimum values of the energy density for a damping loss factors of η = 0.01 and 0.06 are 61.4/50.1 dB and 67.7/32.6 dB, respectively. The difference in the energy density between the maximum and minimum values is known to increase as the damping loss factor increases, which suggests that the global variation of the energy density increases with a high damping loss factor.

Figure 7.

Energy density distribution of the elastic cube on quarter of mid-cross section when $\mathrm{f}=700\mathrm{Hz}$ and $\mathrm{\eta }=0.01$ . The reference energy density is $1\times {10}^{-12}\mathrm{J}/{\mathrm{m}}^3$ . ${\mathrm{L}}_{\mathrm{x}}={\mathrm{L}}_{\mathrm{y}}={\mathrm{L}}_{\mathrm{z}}=200\mathrm{m}$ .

Figure 8.

Energy density distribution of the elastic cube on quarter of mid-cross section when $\mathrm{f}=700\mathrm{Hz}$ and $\mathrm{\eta }=0.06$ . The reference energy density is $1\times {10}^{-12}\mathrm{J}/{\mathrm{m}}^3$ . ${\mathrm{L}}_{\mathrm{x}}={\mathrm{L}}_{\mathrm{y}}={\mathrm{L}}_{\mathrm{z}}=200\mathrm{m}$ .

Conclusion

This study newly developed the vibrational energy governing equations for elastic 3D solid in order to extend the range of applications of the EFA to 3D vibrational problems. The energy density and intensity in 3D elastic solids, obtained from the exact elasticity equation for the dilatational wave in terms of the scalar potential function, is calculated to derive the energy differential equation for the dilatational wave. Time- and space-averaging were applied to derive the energy density and intensity assuming that the structural damping is small. The energy governing equation for a 3D elastic solid is derived by developing the relationship between the time- and space-averaged energy density and intensity.

The developed energy flow model for dilatational wave was applied to the analysis of vibrating cube-shaped structure by a spherical cavity source. And the developed energy governing equation was verified by comparing its results to those obtained with NASTRAN. A comparison of the results for the energy density distribution shows that the approximate energy density calculated using the derived energy flow model represents the global variation of the vibrational response quite well, with reliable results. In addition, as the exciting frequency and structural damping increase, the developed energy flow solutions spatially decrease more quickly.

All these taken together, when the developed energy flow model is applied to the 3D structure, it can analyze the vibration problem of the high-frequency range which is impossible with the conventional FEM and boundary element method (BEM). Also, if energy flow model for distortion waves is developed in the future, it will be possible to extend the application range of EFA to arbitrary shape structures having general boundary conditions.