Prediction of the Transient Local Energy by Energy Finite Element Analysis

Abstract

Energy finite element analysis (EFEA) has been successfully applied to steady-state response prediction over the past three decades. Compared with other energy-based methods, such as statistical energy analysis (SEA), EFEA can consider more local structural information without increasing the computational consumption too much, which makes it attractive. Inspired by the transient local energy approach (TLEA), a general transient energy balance equation was derived by assuming that the plane wave condition is satisfied. The properties of the energy balance equation were studied, and the analytical solutions with different initial conditions were provided. Utilizing the derived transient energy balance equation, transient EFEA is proposed, which has the same advantages as EFEA. A general formula is presented for the energy transmission coefficients of any number of coupled in-plane beams. The present approach was validated using a single beam and a coupled collinear beam structure under unloading conditions. The coupled collinear beams were also investigated using constant and quasi-static input power. The validation results show that TEFEA can accurately predict the local response of the structure. All of these results were compared with those of finite element analysis (FEA), simplified TEFEA (sTEFEA), transient statistical energy analysis (TSEA), and analytical formulas.

1. Introduction

The dynamic response prediction of complex structures has always been an important topic in product design. Generally, finite element analysis (FEA) is preferred because it is straightforward, proven, and widely used. However, FEA becomes less accurate and time-consuming as the frequency and structural complexity increase. Additionally, the uncertainties in the material properties and structural connections due to machine or assembly error make FEA modeling more difficult [1].

Considering these limitations of FEA, some alternative methods have been proposed [2], of which improved deterministic and energy-based methods are well-known. There are many improvements to deterministic methods. For example, the meshless method [3,4] aims to make the model independent of the mesh, which reduces some computation and make it applicable to more complex structures. The dynamic stiffness method (DSM) [5] can always provide more accurate results independent of the number of elements. Recently, some new approaches, such as deep learning and neural network techniques [6,7], have also attracted widespread attention owing to their high efficiency and advantages in addressing the uncertainties. Compared to deterministic methods, energy-based methods are more efficient and statistical, among which statistical energy analysis (SEA) is probably the most popular. The subsystem is the basic component of SEA, such that SEA models are usually simpler than FEA models. Some benefits of statistical properties and energy averaging are that uncertainties can be avoided, even with limited knowledge of the system details. However, SEA also has its flaws. Firstly, its results are valid only if the number of modes is sufficiently large; thus, SEA is usually more appropriate for high-frequency systems. Secondly, all of the SEA parameters are averaged over the entire subsystem. The variations in the material parameters or outer excitations in the subsystem cannot be considered, making the results less accurate. Furthermore, SEA cannot provide the energy distribution inside the subsystem [8].

Transient SEA (TSEA) has a long history that can be traced back to the formulas of Manning and Lee [9]. Furthermore, Pinnington and Lednik considered systems with two oscillators [10] and coupled beams [11], both excited with an impulse. Recently, Langley [12,13] derived the TSEA equations by employing the Priestley description of a non-stationary vibration, and applied them to coupled plates with an impulse. However, similarly to SEA, TSEA can only obtain the energy of the entire subsystem, and will lose the detailed energy distribution information.

EFEA is an alternative to SEA, and is applicable to high-frequency vibration analysis. As an energy-based method, it also requires averaging the energy density over space and frequency. However, unlike SEA, the energy density in EFEA is averaged over one wavelength, so the locally averaged energy distribution can still be obtained. The local effects, such as non-uniform damping and spatially varying loads or geometries, can also be considered compared to SEA. Therefore, EFEA has attracted the attention of researchers since its introduction. Originally, its energy density governing equation was established via an analogy with heat conduction [14]. Based on the assumption of uncorrelated plane waves, EFEA was quickly extended to rods [15], beams [15], members [16], and plates [17,18,19]. Some complex theories, such as the Timoshenko beam [20], the Mindlin plate [21,22], and the Rayleigh–Love and Rayleigh–Bishop rod theories [23], have also been adopted by EFEA. It has also been applied to some built-up structures, such as composite structures [24], beams with stepped thickness and variable cross-sections [25], and stiffened plates [26,27,28,29]. An energy flow analysis model for free-layer damped treated panels has recently been developed [30]. External factors, such as thermal loads [31] and external forces [32], have also been considered. EFEA has also been validated by applying it to some complex structures, such as aircraft [33,34,35], ships [27,36], and vehicles [37,38,39]. Although EFEA has been successfully used, few studies have been conducted on transient analysis based on EFEA. Ichchou developed a new transient equation called the transient local energy approach (TLEA) [40]. Compared with TSEA, TLEA can achieve better rise time and peak energy level results when applied to a two-oscillator system [41]. TLEA has also been used to predict the shock response of a cantilever beam [42]. Chen applied TLEA and TSEA to complex systems and found that TLEA had the advantages of computational efficiency over FEA, and was more accurate than TSEA [43].

In previous studies [41,42], TLEA was only applied to simple structures, and only limited initial conditions were studied. In addition, using steady-state energy intensities in boundary conditions is often incorrect for transient systems. EFEA is mainly used for steady-state response analyses. Therefore, transient EFEA (TEFEA) is proposed based on the transient local equation in this study. It can be used for more complex structures and shares the same advantages as EFEA, such as being more computationally efficient than FEA and containing more structural local information than TSEA. Structures with many uncertainties (such as boundary conditions, material properties, and manufacturing errors) can benefit from the use of this approach. In Section 2.1, inspired by TLEA, a more general method is utilized to derive the transient energy flow equation, and the basic assumptions of TEFEA are discussed. Section 2.2 describes analytical solutions to the transient local energy equation for one-dimensional structures. Accordingly, two common initial conditions are introduced. In Section 3, the newly developed approach is validated through comparison with FEA, TSEA, and a simplified TEFEA (sTEFEA) using diffuse field approximation. In Appendix A and Appendix B, a method similar to EFEA is used to derive the energy flow relationship and the finite element equation.

2. Energy Flow Model

This section uses a more general method to derive the transient energy flow equation that can be used in multidimensional structures. Some basic assumptions that should be satisfied are the following: (i) the system can be simplified as a hysteretic damping system, and the damping force is weak; (ii) evanescent waves (near-field displacement) are omitted; (iii) the interference between different traveling waves is ignored; (iv) the waves in the structure can be simplified as plane waves; and (v) the average energy and average power input vary significantly slower than the oscillations of the system [12].

2.1. Energy Density Governing Equation

First of all, the energy balance (or power flow balance) relationship should be introduced, as follows [14]:

(1)$\frac{\partial E}{\partial t}+\nabla \cdot \overrightarrow{I}+{\pi }_{diss}={\pi }_{in},$

where $E$ is the energy density, $\overrightarrow{I}$ is the energy intensity, ${\pi }_{diss}$ is the energy dissipated by the system damping, and ${\pi }_{in}$ is the input power. If assumptions (i) and (iii) are satisfied, the energy dissipation can be expressed as follows [44]:

(2)${\pi }_{diss}=\eta \omega E,$

where $\eta$ is the structural damping coefficient and $\omega$ is the circle frequency.

The equation based on the diffuse field can easily be derived from the steady-state energy flow equation [15]. At time $t$ and position $r$, the displacement of a unidirectional plane traveling wave can be written as follows:

(3)$w\left(r,t\right)=A\cdot e^{-K\cdot r}{e}^{i\omega t},$

where $K$ is the wave number. Taking flexure waves as an example, the energy density and energy intensity can be expressed in terms of the displacement as:

(4)$\begin{array}{l}E=\frac{\rho }{2}\frac{\partial w}{\partial t}\cdot {\left(\frac{\partial w}{\partial t}\right)}^{*}=\frac{D}{2}\frac{{\partial }^{2}w}{\partial x^2}\cdot {\left(\frac{{\partial }^{2}w}{\partial x^2}\right)}^{*}=\frac{\rho {\omega }^2}{2}{\left|A\right|}^2{e}^{-\frac{\eta k_1}{2}\cdot x}=\frac{D{\left|K\right|}^4}{2}{\left|A\right|}^2{e}^{-\frac{\eta k_1}{2}\cdot x}\\ I=-\frac{D}{2}\left[\frac{{\partial }^{2}w}{\partial x^2}\cdot {\left(\frac{{\partial }^{2}w}{\partial t\partial x}\right)}^{*}-\frac{{\partial }^{3}w}{\partial x^3}\cdot {\left(\frac{\partial w}{\partial t}\right)}^{*}\right]=D\omega K^2{k}_1{\left|A\right|}^2{e}^{-\frac{\eta k_1}{2}\cdot x}\end{array},$

where $D=Y\cdot h^3/12$ is the stiffness parameter, $Y$ is Young’s modulus, $h$ is the thickness of the structure, and $k_1$ is the real part of the wavenumber. If the damping is weak, the relationship between the energy density and intensity can easily be derived as follows:

(5)$I=-\frac{C_g{}^2}{\eta \omega }\frac{\partial E}{\partial x}.$

Substituting Equation (5) into Equation (1), one can yield:

(6)$\frac{\partial E}{\partial t}-\frac{C_g^2}{\eta \omega }\frac{{\partial }^{2}E}{\partial x^2}+\eta \omega E={\pi }_{in}.$

If reflected waves exist and assumption (ii) is satisfied, the displacement of the plane wave should be the following:

(7)$w\left(r,t\right)=\left(A\cdot e^{-K\cdot r}+B\cdot e^{K\cdot r}\right)e^{i\omega t}.$

Therefore, the energy density and energy intensity should be as follows:

(8)$\begin{array}{l}E=\frac{1}{4}\left(D{\left|K\right|}^4+\rho \cdot {\omega }^2\right)\left({\left|A\right|}^2\cdot e^{-\frac{\eta k_1}{2}x}+{\left|B\right|}^2\cdot e^{\frac{\eta k_1}{2}x}+2\mathrm{Re}\left(A{B}^{*}\cdot e^{-i\cdot 2k_1\cdot x}\right)\right)\\ I=D\omega K^2{k}_1\left({\left|A\right|}^2\cdot e^{-\frac{\eta k_1}{2}x}-{\left|B\right|}^2\cdot e^{\frac{\eta k_1}{2}x}\right)+\frac{1}{2}D\eta \omega K^2{k}_1\cdot \mathrm{Im}\left(A{B}^{*}\cdot e^{-i\cdot 2k_1\cdot x}\right)\end{array}$

where $\mathrm{Re}\left(\ast \right)$ and $\mathrm{Im}\left(\ast \right)$ in Equation (8) are the fluctuation items, and are eliminated if averaged over one wavelength. Using Equation (8), a similar equation to Equation (5) can also be derived as follows:

(9)$⟨I⟩=-\frac{C_g{}^2}{\eta \omega }\frac{\partial ⟨E⟩}{\partial x},$

where $⟨\ast ⟩$ indicates that the value is averaged over one wavelength.

According to references [17,18,24], the equation for a two-dimensional plane-wave system, such as a plate, can be obtained as:

(10)$⟨\overrightarrow{I}⟩=-\frac{C_g^2}{\eta \omega }\nabla ⟨E⟩.$

Substituting Equation (10) into Equation (1) can yield the following:

(11)$\frac{\partial ⟨E⟩}{\partial t}-\frac{C_g^2}{\eta \omega }{\nabla }^2⟨E⟩+\eta \omega ⟨E⟩={\pi }_{in}.$

Equations (6) and (11) explicitly require that the system be in a quasi-steady state, such that the energy intensity can be represented by Equations (5) and (9). The quasi-steady state is also required for TSEA [12].

Notably, the energy density in Equations (9)–(11) is averaged over one wavelength. Additionally, for systems of limited size, the energy density and energy intensity in EFEA are also averaged over a specific frequency band [1], owing to the presence of the reflected waves. TEFEA should also meet these requirements. Similar formulas can also be derived for longitudinal and shear waves [15,19].

Another energy intensity and energy density relationship stated in [44] can be extended to multidimensional plane wave systems (Appendix B, the case without reflected waves):

(12)$\overrightarrow{I}=\overrightarrow{C_g}\cdot E.$

Equation (12) is derived based on diffuse wave fields, and cannot be used directly if reflected waves exist. However, if the energy density is decomposed in all directions like the wave number, a similar relationship in Equation (A28) can be obtained. The detailed relationship can be found in Appendix B. Also, the energy intensity and density should be averaged over one wavelength, as in Equations (8)–(11), and the symbol $⟨*⟩$ is omitted below. If there is no external force, Equation (1) can be simplified as:

(13)$\frac{\partial E}{\partial t}+\nabla \cdot \overrightarrow{I}+\eta \omega E=0.$

For two-dimensional structures, it can be rewritten as:

(14)$\frac{\partial E}{\partial t}+{\displaystyle \sum _i\frac{\partial I_i}{\partial x_i}}+\eta \omega E=0,\left(i=x,y\right).$

Substituting Equations (A26) and (A29) into Equation (14), the energy flow balance equation in every direction can be expressed as:

(15)$\frac{\partial E_i^{\pm }}{\partial t}+\frac{\partial I_i^{\pm }}{\partial y}+\eta \omega E_y^{\pm }=0,\left(i=x,y\right),$

where the superscript “+” indicates the energy of the traveling wave along the positive direction of the axis, and superscript “–“ indicates the energy in the opposite direction.

The energy flow balance relationship in the positive direction of the y-axis can be written as:

(16)$\frac{\partial E_y^{+}}{\partial t}+\frac{\partial I_y^{+}}{\partial y}+\eta \omega E_y^{+}=0.$

Substituting Equation (A28) into Equation (16) can yield:

(17)$\frac{\sin \theta }{C_g}\frac{\partial I_y^{+}}{\partial t}+\frac{C_g}{\sin \theta }\frac{\partial E_y^{+}}{\partial y}+\eta \omega \frac{\sin \theta }{C_g}{I}_y^{+}=0.$

Equation (16) is differentiated with respect to time t, Equation (17) is differentiated with respect to y, and then the items of energy intensity $I^{+}$ are eliminated using Equation (17). The energy density governing equation can be written as:

(18)$\frac{{\partial }^2{E}_y^{+}}{\partial t^2}+2\eta \omega \frac{\partial E_y^{+}}{\partial t}+{\left(\eta \omega \right)}^2{E}_y^{+}-\frac{C_g^2}{{\sin }^2\theta }\frac{{\partial }^2{E}_y^{+}}{\partial y^2}=0.$

Summing the energy in both directions and applying Equation (A25), the energy density balance equation for the y-axis can be obtained as:

(19)$\frac{{\partial }^2{E}_y}{\partial t^2}+2\eta \omega \frac{\partial E_y}{\partial t}+{\left(\eta \omega \right)}^2{E}_y-C_g^2\frac{{\partial }^{2}E}{\partial y^2}=0.$

The same energy relationship exists on the x-axis, and the total energy density balance equation can be expressed as:

(20)$\frac{{\partial }^{2}E}{\partial t^2}+2\eta \omega \frac{\partial E}{\partial t}+{\left(\eta \omega \right)}^{2}E-C_g^2\cdot {\nabla }^{2}E=0.$

If the external power input ${\pi }_{in}$ exists, the governing equation can be rewritten as:

(21)$\frac{1}{\eta \omega }\frac{{\partial }^{2}E}{\partial t^2}+2\frac{\partial E}{\partial t}+\eta \omega E-\frac{C_g^2}{\eta \omega }{\nabla }^{2}E={\pi }_{in}.$

Equation (21) has the same form as that of TLEA in [40]. Comparing Equations (1) and (21), the divergence of the energy intensity can be obtained as:

(22)$\nabla \cdot \overrightarrow{I}=\frac{1}{\eta \omega }\frac{{\partial }^{2}E}{\partial t^2}+\frac{\partial E}{\partial t}-\frac{C_g^2}{\eta \omega }{\nabla }^{2}E.$

For example, if only the steady state is considered, the time-dependent items in Equations (21) and (22) can be neglected. Thus, the relationship between the energy intensity and energy density, found in [15,16,17,18,19], can be expressed as the following Equations (23) and (24):

(23)$\nabla \cdot \overrightarrow{I}=-\frac{C_g^2}{\eta \omega }{\nabla }^{2}E$

(24)$-\frac{C_g^2}{\eta \omega }{\nabla }^{2}E+\eta \omega E={\pi }_{in}$

Like EFEA, TEFEA equations can be established using finite element techniques (Appendix A). Equations (6) and (11) are TSEA-like equations that can be formalized as Equation (A11). This simplified equation is abbreviated as sTEFEA in context.

Although Equations (13)–(21) assume the existence of reflected waves, they are also applicable to structures where direct fields dominate (corresponding to Equation (A12)). From another perspective, a multidimensional plane wave can be viewed as a plane wave with amplitude modulation in other dimensions. The energy density and energy intensity of a multidimensional plane wave are the superpositions of their components [17,18,19,24].

In addition, assumption (iii) enables the energy to be expressed as the superposition of the energy per wave, and assumption (v) enables the input power to be described by the amplitude of the external force and the admittance of the input point. Also, Equations (11) and (21) require the energy density to be averaged over one wavelength if reflected waves exist. As mentioned in [1], if the structure size is limited, the results of this transient method should be averaged over a certain frequency range. These requirements are identical to those of EFEA, and are the basis for applying the transient local energy equation to TEFEA.

2.2. Analytical Solution of the Transient Equation

If only one-dimensional plane waves are considered, such as Euler–Bernoulli beams, Equation (20) can be simplified as:

(25)$\frac{{\partial }^{2}E}{\partial t^2}+2\eta \omega \frac{\partial E}{\partial t}+{\left(\eta \omega \right)}^{2}E-C_g^2\cdot \frac{{\partial }^{2}E}{\partial s^2}=0.$

Equation (26) can be solved by separating the variables. Suppose that the general solution takes the following form:

(26)$E\left(s,t\right)=f\left(t\right)g\left(s\right).$

Substituting it into Equation (25) yields:

(27)$f^{″}+2\eta \omega f^{\prime }+\left(1-m^2\right){\left(\eta \omega \right)}^{2}f=0$

(28)${\left(\frac{C_g}{\eta \omega }\right)}^2\frac{{\partial }^{2}g}{\partial s^2}-m^{2}g=0$

where $m$ is a parameter related to the spatial distribution of the energy density, and its value is related to the initial conditions of Equation (25).

The general solutions to Equations (27) and (28) can be expressed as:

(29)$f\left(t\right)=e^{-\eta \omega t}\left(A_1{e}^{-m\eta \omega t}+A_2{e}^{m\eta \omega t}\right)$

(30)$g\left(s\right)=B_1{e}^{-\frac{m\eta \omega }{C_g}\cdot s}+B_2{e}^{\frac{m\eta \omega }{C_g}\cdot s}.$

From Equation (22), it can be seen that the energy intensity is also a time-dependent value. However, in some special cases, the energy intensity can be simplified. For example, when a pulse impinges on a beam [42], stationary boundary conditions can be applied unless the pulse reaches system boundaries. Assuming that the length of the beam is $l$, the boundary condition can be expressed as:

(31)$\begin{array}{l}I\left(t,0\right)={\left.-\frac{C_g^2}{\eta \omega }\frac{\partial e}{\partial s}\right|}_{s=0}=0,\\ I\left(t,l\right)={\left.-\frac{C_g^2}{\eta \omega }\frac{\partial e}{\partial s}\right|}_{s=l}=0.\end{array}$

If $m^2<0$, Equations (29) and (30) are applied to Equation (31), and the unknown parameters can be obtained as:

(32)$\begin{array}{l}{B}_1=B_2\\ m=\pm \cdot i\cdot \frac{C_g}{\eta \omega }\frac{n\pi }{l}\end{array}.$

Thus, the energy density can be expressed as:

(33)$E\left(s,t\right)=C_0+e^{-\eta \omega t}{\displaystyle \sum _n{C}_n\cos \left(\frac{n\pi C_g}{l}t\right)\cos \left(\frac{n\pi }{l}s\right)}.$

The value $C_n\left(n=0\dots \infty \right)$ depends on the input power. This solution can be found in [42], which describes the energy density history of a cantilever beam subjected to a transverse unit impulse at its free end.

If $m^2=0$, Equation (27) has two duplicate eigenvalues. Considering that the energy density is attenuated with time owing to the non-zero loss factor, only Equation (34) is valid. The spatial distribution of the energy density will be uniform.

(34)$f\left(t\right)=C{e}^{-\eta \omega t}.$

If $m^2>0$, the value of $m$ will be real, and every item in Equation (30) will be exponential. Specifically, if $m^2=1$, the corresponding energy density distribution is the same as that in the steady state, and the solution of Equation (27) can be obtained as:

(35)$f\left(t\right)=C_0+C_1{e}^{-2\eta \omega t}.$

Given the initial conditions, applying Equation (35) to Equation (2) yields:

(36)$\left\{\begin{array}{c}f\left(0\right)=C_0+C_1=1\\ \dot{f}\left(0\right)=-2\eta \omega C_1=-\eta \omega \end{array}\right.\Rightarrow \left\{\begin{array}{c}{C}_0=\frac{1}{2}\\ C_1=\frac{1}{2}\end{array}\right..$

Equation (36) can describe a situation in which a beam initially has a steady load, which is removed after it finally settles into a steady state. However, immediate removal of the load would cause disturbances on the beam. Both FEA and TEFEA can predict this phenomenon, as described in Section 3.

From the previous analysis, it can be known that Equation (25) has different forms of solutions according to the initial energy distribution. Reference [42] only includes the condition $m^2<0$, and other cases are considered in this study.

2.3. Energy Transmission Coefficients for Beams

According to Appendix A, energy transmission coefficients are required to establish the coupling matrix. These coefficients can be obtained from [45], which provides the energy transmission formulas for some built-up structures (beam-to-beam, plate-to-plate, and structure-to-acoustic field). Furthermore, from [46], a general method is used to obtain the energy transmission coefficients of any number of coupled in-plane beams.

According to the Euler–Bernoulli beam theory, the longitudinal and transverse motion equations of the beam are expressed as:

(37)$\begin{array}{l}Y\cdot A\frac{{\partial }^{4}u}{\partial x^4}=\rho \cdot A\frac{{\partial }^{2}u}{\partial t^2}\\ D\frac{{\partial }^{4}w}{\partial x^4}=-\rho \cdot A\frac{{\partial }^{2}w}{\partial t^2}\end{array},$

where $u$ and $w$ represent the longitudinal and transverse displacements, respectively; $D=Y\cdot b\cdot h^3/12$ is the stiffness parameter; $Y$ is Young’s modulus; $b$ and $h$ are the width and height of the cross-section, respectively; $\rho$ is the mass per unit volume; and $A$ is the cross-sectional area.

Suppose there are $i$ beams that are linked at a joint (Figure 1). The force acting on the cross-section can be written in the form of Equation (38):

(38)$N=Y\cdot A\frac{\partial u}{\partial x},S=-D\frac{{\partial }^{3}w}{\partial x^3},M=D\frac{{\partial }^{2}w}{\partial x^2}.$

By analyzing the forces exerted on the joint, Equation (39) can be readily obtained:

(39)$\left\{\begin{array}{l}{\displaystyle \sum _i{N}_i\cos {\theta }_i-S_i\sin {\theta }_i}=0\\ {\displaystyle \sum _i{N}_i\sin {\theta }_i+S_i\cos {\theta }_i}=0\\ {\displaystyle \sum _i{M}_i}=0\end{array}\right..$

The displacement of the joint can be represented by the local displacement of beam i at $\mathrm{x}=0$ as follows:

(40)$\left\{\begin{array}{l}u=u_i\cos {\theta }_i-w_i\sin {\theta }_i\\ w=u_i\sin {\theta }_i+w_i\cos {\theta }_i\\ \partial w/\partial x={\partial w_i/\partial x}_i\end{array}\right..$

The force and displacement in the local coordinates at the joint of beam i can be expressed as:

(41)$Q_i={\left\{N_i,S_i,M_i\right\}}^T,U_i={\left\{u_i,w_i,\partial w_i/\partial x_i\right\}}^T.$

Therefore, Equations (39) and (40) can be rearranged as Equation (43), where $R_i$ is the coordinate transformation matrix, and $U_0$ represents the displacement of the joint at the global coordinate:

(42)$R_i=\left[\begin{array}{ccc}\cos {\theta }_i& -\sin {\theta }_i& 0\\ \sin {\theta }_i& \cos {\theta }_i& 0\\ 0& 0& 1\end{array}\right]$

(43)${\displaystyle \sum _i{R}_i{Q}_i}=0,U_0=R_i{U}_i$

The longitudinal and transverse displacements of beam $i$ are assumed to be of the following form:

(44)$\begin{array}{l}{u}_i={\alpha }_{Li}{e}^{{\mu }_{iL}{x}_i}\\ w_i={\alpha }_{F1i}{e}^{{\mu }_{F1i}{x}_i}+{\alpha }_{F2i}{e}^{{\mu }_{F2i}{x}_i}\end{array},$

where subscripts $i$, $L$, and $F$ denote beam $i$, the longitudinal wave, and the flexural wave, respectively. ${\mu }_{iL}$, ${\mu }_{F1i}$, and ${\mu }_{F2i}$ are the longitudinal wave number and two flexural wave numbers, respectively. The wave number of the reflected and transmitted waves can be obtained as follows:

(45)$\left\{\begin{array}{l}{\mu }_{Li}=-k_{Li}\\ {\mu }_{F1i}=-k_{Fi}\\ {\mu }_{F2i}=-j{k}_{Fi}\end{array}\right.,k_{Li}=\omega \sqrt{\frac{{\rho }_i}{Y_i}},k_{Fi}=\sqrt[4]{\frac{{\rho }_i{A}_i{\omega }^2}{D_i}}.$

Substituting Equations (38) and (44) into Equation (41), the force $Q_i$ can be rewritten as:

(46)$Q_i=\left[\begin{array}{ccc}{\left(EA\right)}_i{\mu }_L& 0& 0\\ 0& -D_i{\mu }_{F1i}^3& -D_i{\mu }_{F2i}^3\\ 0& D_i{\mu }_{F1i}^2& D_i{\mu }_{F2i}^2\end{array}\right]\left[\begin{array}{ccc}1& 0& 0\\ 0& \frac{{\mu }_{F2i}}{{\mu }_{F2i}-{\mu }_{F1i}}& \frac{-1}{{\mu }_{F2i}-{\mu }_{F1i}}\\ 0& \frac{-{\mu }_{iF1}}{{\mu }_{F2i}-{\mu }_{F1i}}& \frac{1}{{\mu }_{F2i}-{\mu }_{F1i}}\end{array}\right]U_i.$

For simplicity, Equation (46) can be written as:

(47)$Q_i=T_i\cdot U_i,$

where

(48)$T_i=\left[\begin{array}{ccc}{T}_{i,11}& 0& 0\\ 0& T_{i,22}& T_{i,23}\\ 0& T_{i,32}& T_{i,33}\end{array}\right]$

and

(49)$\begin{array}{l}{T}_{i,11}=E_i{A}_i{\mu }_{Li},\hspace{1em}{T}_{i,22}=D_i{\mu }_{F1i}{\mu }_{F2i}\left({\mu }_{F1i}+{\mu }_{F2i}\right)\\ T_{i,23}=-D_i\left({\mu }_{F2i}^2+{\mu }_{F1i}{\mu }_{F2i}+{\mu }_{F1i}^2\right)\\ T_{i,32}=-D_i{\mu }_{F1i}{\mu }_{F2i},\hspace{1em}\hspace{1em}{T}_{i,33}=D_i\left({\mu }_{F2i}+{\mu }_{F1i}\right)\end{array}.$

At the joint, the displacement can be expressed as:

(50)$U_i=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 1\\ 0& {\mu }_{F1}& {\mu }_{F2}\end{array}\right]\left\{\begin{array}{c}{\alpha }_{iL}\\ {\alpha }_{iF1}\\ {\alpha }_{iF2}\end{array}\right\}.$

Equation (47) is valid only if the incident wave is not on plate $i$. Otherwise, it is the result of adding the reflected waves to the incident waves. Thus, if $Q_m$ and $U_m$ denote the force and displacement superposition of the incident and reflected waves, Equation (47) can be expressed as:

(51)$Q_m-Q_m^{\prime }=T_m\cdot \left(U_m-U_m^{\prime }\right),$

where $U_m^{\prime }$ and $Q_m^{\prime }$ represent the incident wave displacement and force, respectively. Equation (51) can also be rewritten as:

(52)$\begin{array}{l}{Q}_m=T_m\cdot U_m-Q_{m,inc}\\ Q_{m,inc}=T_m\cdot U_m^{\prime }-Q_m^{\prime }\end{array}.$

Substituting Equation (52) into Equation (43), the balanced relationship at the joint can be expressed as:

(53)${\displaystyle \sum _i{R}_i\cdot T_i\cdot R_i^T{U}_0}=R_m{Q}_{m,inc}.$

Assuming that a flexural wave of unit amplitude on beam m is incident to the joint, the reflected and transmitted waves will propagate along the positive direction of the local coordinate $x_i$; thus, the flexural and longitudinal wave numbers must be negative (including the negative imaginary value). However, the incident wave has the opposite direction, so positive wavenumbers are required. From Equations (38) and (44), the force and displacement due to the incident wave can be expressed as:

(54)$Q_m^{\prime }={\left\{0,j{E}_m{I}_m{k}_{Bm}^3,-E_m{I}_m{k}_{Fm}^2\right\}}^T,U_m^{\prime }=\{0,1,j{k}_{Fm}\}.$

Substituting Equation (54) into Equation (53), the displacement $U_0$ can be obtained. The local displacement of beam $i$ can be obtained by solving Equation (43). The amplitude of each wave can be obtained by solving Equation (50).

The energy intensity of the longitudinal and flexural wave is as follows:

(55)$P_L=\frac{1}{2}{\alpha }_L^2\rho {\omega }^3/k_L,P_F={\alpha }_F^2\rho {\omega }^3/k_F.$

The transmission coefficient of the wave of type S from beam $i$ to the wave of type $T$ from beam $j$ can be written as:

(56)${\tau }_{STij}=\frac{P_{Tj}}{P_{Si}}.$

3. Validation and Discussion

The present approach uses the basic assumptions of EFEA, so it also inherits the limitations of EFEA. In addition, the coupling loss factor (CLF) in Section 2.3 is derived in a steady state, and sometimes it will be inaccurate [47]. However, if the energy does not change much over a short time, a steady-state CLF can also be used. The Euler–Bernoulli beam model is adopted because it represents a single bending wave. TEFEA using other beam theories can be extended by combining different types of waves [20].

In this section, the derived TEFEA equations are validated numerically by applying them to the beam system. The simulation results of TEFEA are compared with those of FEA, sTEFEA, and TSEA.

Two initial conditions are considered: unloading and loading. In this context, the unloading condition refers to the immediate load removal after the system reaches a steady state. In particular, a steady point force is initially applied at the specified location on the beam, and is removed after the system has been in a steady state for a long enough time. Equation (24) is used to calculate the steady-state energy density distribution of EFEA. The steady-state result is applied to the system as the initial condition of the transient analysis. The loading condition indicates that a continuous force is exerted on the beams with zero initial energy. Two types of loads are investigated in this section.

sTEFEA has the same initial conditions as TEFEA. Equation (35) gives the analytical results with the initial energy of FEA.

All of the transient results of FEA are averaged over one wavelength. The Newmark–Beta method is used to solve the TEFEA equations. The sTEFEA and TSEA equations are solved using the fourth-order Runge–Kutta formulas.

3.1. Numerical Simulations and Verifications

The response of a single beam pinned at both ends was investigated to check the convergence of Equations (11) and (21). Its material parameters are as follows: density $2700 \mathrm{k}\mathrm{g}/{\mathrm{m}}^3$, elastic modulus ${71\times 10}^9 \mathrm{N}/{\mathrm{m}}^2$, and Poisson’s ratio $0.33$. The geometric properties of the beam are as follows: cross-section $b\times h=0.02 \mathrm{m} \times 0.002 \mathrm{m}$ and entire length $L=5 \mathrm{m}$. The damping factor is 0.01. A harmonic point force with an amplitude of 10 N and a frequency of 4000 Hz is applied at the center of the beam. Figure 2 lists some response points that will be used later. It should be noted that TEFEA, like EFEA, is not sensitive to boundary conditions [1]; thus, only pinned–pined boundary conditions are used.

In Figure 3, the beam is under the unloading condition and meshed with 8, 40, 200, and 1000 elements separately. The results show that the more elements used, the smoother the response curve. An approximate result can also be obtained with 8 elements for TEFEA. The error can be almost ignored if the number of elements reaches 200. For sTEFEA, different mesh sizes have little effect on the simulation results.

A similar verification was performed by adding a varying load (as in Figure 4), and the simulation results are shown in Figure 5. All of the TEFEA results in Figure 3 and Figure 5 were obtained by direct numerical integration.

Equation (A10) can also be solved using mode decomposition instead of direct integration. In this case, the beam is meshed with 200 elements. The simulation results with 20, 40, 80, and 160 modes were compared with the simulation results of the direct integration method. According to Figure 6 and Figure 7, acceptable results can be obtained when the mode number reaches 80. Unlike FEA, the modes in EFEA cannot represent some form of energy vibration, since the energy should not be negative. Therefore, the wavelength criterion in FEA cannot also be used. The time history of the energy in EFEA is related to the envelope of the velocity time history in FEA. Consequently, the energy varies very slowly compared to the displacement or the velocity used in FEA, and an approximate energy result can be obtained using only a limited number of modes. In Section 3.2, 200 elements and 160 modes are used.

3.2. A Single Beam under the Unloading Condition

In this section, the same beam depicted in Figure 2 was used to validate Equations (11) and (21) under the unloading condition. The FEA model was meshed densely such that the modes under 10,000 Hz were valid. Only the transverse modes were considered, and an ensemble of 315 modes below 10,000 Hz was used to calculate the steady-state response. Like EFEA, the steady-state result was used as the initial conditions for FEA. The FEA model was meshed with 2000 elements to ensure accuracy.

Figure 8 shows the steady-state response results for all of the methods, revealing that the EFEA results agree well with the FEA results. However, SEA can only provide an average energy. Therefore, the energy densities of the SEA at the center and two ends of the beam are significantly different from those of FEA, resulting in inaccurate initial energy when TSEA is initialized.

Figure 9 shows the energy density time histories of all of the selected positions. An obvious correlation between FEA and TEFEA can be found, especially at locations far from the excitation source. The dotted lines in Figure 9 indicate the time required for the energy to travel from the excitation source to the corresponding location. Significant changes in the FEA and TEFEA results can be found at the times indicated by the dotted lines. These changes are caused by removing the external force, so the system suddenly becomes unsteady. For TEFEA, the immediate removal of the external input power is equivalent to applying a transient input power of equal magnitude and opposite direction at the location of the force. This transient unsteady state introduces a disturbance propagating to the entire structure at wave group speed. The propagation of this perturbation energy is only predicted by FEA and TEFEA. Compared with FEA, the analytical results (Equation (35)) in Section 2.3 provide a perfect approximation before falling into an unsteady state. This phenomenon is consistent with that reported in [42].

A straightforward comparison can be found in Figure 10, from which changes in the energy spatial distribution can be observed. The FEA and TEFEA results are in good agreement at all of the specified times. Figure 10 also shows that the energy density distribution tends to be uniform over time, corresponding to $m^2=0$ for TEFEA described in Section 2.3.

Figure 9 and Figure 10 also demonstrate that sTEFEA is similar to TSEA. Although they have different initial energy densities, they become consistent over time. For sTEFEA, the beam energy distribution at every location changes immediately after the disturbance, indicating that the disturbance propagates at an infinite speed [42]. This is impossible for vibrational energy propagation in structures. The difference between TEFEA and sTEFEA can be explained using Equations (16) and (17). By adding Equation (17), the relationship between energy density and intensity (Equations (12) and (B17)) is always guaranteed.

The simulations were run on a personal computer with an Intel core i7 4770 processor and 64 GB RAM. The consumption times obtained for FEA, TEFEA, and SEA are 1.50 s, 0.0705 s, and 0.035 s, respectively (average values of 10 runs). All of these methods use the same time step. TSEA is the most efficient, but cannot correctly predict the local energy. Compared to FEA, TEFEA is more efficient and sufficiently accurate.

3.3. Coupled Beams under Unloading Condition

In this section, the coupled system illustrated in Figure 11 is investigated to validate TEFEA and sTEFEA. The two beams are collinear, connected directly at the coupled ends, and pinned at the other ends. Both beams share the same material parameters: density $2700 \mathrm{k}\mathrm{g}/{\mathrm{m}}^3$, elastic modulus ${71\times 10}^9 \mathrm{N}/{\mathrm{m}}^2$, and Poisson’s ratio $0.33$. The size of the left beam is as follows: cross-section $b\times h=0.05 \mathrm{m} \times 0.004 \mathrm{m}$ and length $L=5 \mathrm{m}$. The size of the right beam is as follows: cross-section $b\times h=0.08 \mathrm{m} \times 0.006 \mathrm{m}$ and length $L=5 \mathrm{m}$. Both damping factors were set to 0.01.

A random point force with a frequency range of 3500–4500 Hz was applied at the center of the left beam. The amplitude of the force was 10 N. The energy transmission coefficients were calculated using Equations (53)–(56). The FEA results were averaged over 3500–4500 Hz, and the results at the central frequency were used in TEFEA, sTEFEA, and TSEA. Only the transverse modes were considered, and an ensemble of 315 modes under 10,000 Hz was used to calculate the steady-state response. An ensemble of 500 realizations over the frequency range of 3500–4500 Hz was used. There were 23 modes in the frequency range of 3500–4500 Hz.

The steady-state response results shown in Figure 12 were used as the initial energy of the transient analysis. The EFEA results match the FEA results closely, except for the locations of the excitation point and boundaries. Conclusively, the energy transmission coefficients calculated using Equation (53) agree with the FEA. However, the SEA results are unreliable compared to the FEA results, especially for the right beam.

Similar conclusions can be drawn from Figure 13, as described in Section 3.2. The TEFEA results agree well with the FEA results, and they are all consistent with the analytical solution before being affected by the disturbance. FEA and TEFEA can predict the disturbance propagation, whereas sTEFEA and TSEA cannot. In particular, Figure 13e shows the response at the center of the right beam, and the sTEFEA and TSEA results are significantly different from the FEA results. The time histories of the local energy variation for sTEFEA and TSEA are also not consistent with that for FEA.

3.4. Coupled Beams under the Loading Conditions

This section investigates the TEFEA formulas with continuous random excitations. The coupled beams described in Section 3.3 were adopted. The amplitude of the excitation was 10 N, and the system was initially at rest. It was assumed that the load only acts at 3500–4500 Hz; the energy densities at different locations are shown in Figure 14.

From Figure 14 and Figure 15, it can be observed that all of the responses of FEA, TEFEA, sTEFEA, and TSEA are asymptotic towards their steady-state values. Similarly to the results in the steady-state cases, the TEFEA and sTEFEA results provide good predictions of the energy density when the system reaches the steady state, whereas the TSEA results cannot. The results of FEA and TEFEA show that the response at the unforced positions is not triggered as soon as the excitation is exerted, and the energy propagates away from the excitation at the group speed. In contrast, this process is not predicted by sTEFEA or TSEA. A similar phenomenon can also be observed in Figure 15b.

It is also noticeable from Figure 14c that the position at which the force is exerted acquires an instantaneous velocity in an extremely short time, like a shock, because when a high-frequency force is applied to a structure at rest, the structure will behave as if a shock is exerted. In contrast, the TEFEA results change more slowly. TEFEA can only predict the arrival time of the shock, but not its amplitude. sTEFEA and TSEA cannot predict this shock. Alternatively, a quasi-static load (Figure 4) was used to avoid the shock effect. The maximum amplitude of the excitation was 10 N.

If using the load in Figure 4, the initial velocity resulting from the shock effect can be avoided. The results are depicted in Figure 16 and Figure 17, from which good agreement can be found between the FEA, TEFEA, and sTEFEA results, except that the results of TEFEA lag slightly behind those of FEA. As the energy in the structure increases slowly from zero, the structure also has enough time to reach a steady state. No apparent disturbance propagation can be observed in any time history. The FEA, TEFA, and sTEFEA results are in good agreement for such quasi-static loads.

4. Conclusions

This study proposed a transient energy balance equation using a general method, indicating that the transient energy equation could be applied to multidimensional plane wave systems. Similarly to EFEA, a transient form of EFEA was established. The basic assumptions of TEFEA were discussed in Section 2.1. Similarly to steady-state EFEA, TEFEA requires averaging over a wavelength and a sufficient frequency range. Studying the analytical solutions of the transient local energy equation shows that this equation has different solution forms that depend on the initial energy distribution. Based on the analytical solutions, two initial conditions were considered in Section 3. The transmission coefficients of the coupled beams were calculated using the derived formula in Section 2.3.

The TEFEA results were validated using a single beam and coupled collinear beams. These results were compared with those of FEA, TSEA, and sTEFEA under loading and unloading conditions. The transient energy propagation was well predicted by TEFEA under both initial conditions. TEFEA requires a lower mesh density and fewer modes without loss of accuracy for high-frequency vibration analysis, which makes it more computationally efficient than FEA. Compared with sTEFEA and TSEA, TEFEA can predict both the time history of the energy propagation and the space distribution of the energy density under unloading and loading conditions. sTEFEA and TSEA failed to capture the energy propagation process, and the energy distribution changed instantly at all locations, even if the energy disturbance was exerted at a point. However, if the external input power is quasi-static, the energy propagation process is negligible, and sTEFEA can also be used to predict the time history of the diffuse energy.

Although the proposed method can describe the energy time history well, it also has certain limitations. For example, it relies on the same assumptions as the steady-state EFEA, and requires the energy results to be averaged over a specified frequency band [1] and one wavelength (see Section 2.1). Therefore, it is suitable for high-frequency vibration analysis. Additionally, it is assumed that the waves in the structure are plane-wave, which can be satisfied if the reverberant field dominates [48]. Otherwise, the contribution of the direct field should also be included [8].

Furthermore, the energy governing equation can be applied to more complex structures and initial conditions if the initial energy density distribution can be obtained correctly. In addition, TEFEA can also be used for time-varying systems, but the transient energy balance equation needs to be re-derived to consider the time-varying parameters [49].

Author Contributions

Conceptualization, H.W. and K.Y.; methodology, H.W.; software, H.W.; validation, H.W., K.Y. and R.Z.; formal analysis, H.W., K.Y. and R.Z.; investigation, H.W., K.Y. and R.Z.; data curation, H.W.; writing—original draft preparation, H.W.; writing—review and editing, H.W., K.Y. and R.Z.; project administration, H.W. and K.Y. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The data that support the findings of this study are available from the corresponding author, K. Yu, upon reasonable request.

Acknowledgments

The authors would like to thank the editor and the anonymous reviewers for their valuable suggestions on improving this study.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1. Schematic of coupled in-plane beams.

Figure 2. Single beam pinned at both ends.

Figure 3. Responses of beams with different numbers of elements under the unloading condition at (a) P1 and (b) P5. These results were obtained with 8, 40, 200, and 1000 elements.

Figure 4. Normalized amplitude of the load.

Figure 5. Responses of beams with different numbers of elements under the loading condition at (a) P1 and (b) P5. These results were obtained with 8, 40, 200, and 1000 elements.

Figure 6. Responses of beams at (a) P1 and (b) P5 under the unloading condition. These results were obtained with 20, 40, 80, and 160 modes; “direct” indicates the direct integration solution.

Figure 7. Responses of the beam at (a) P1 and (b) P5 under the loading condition. These results were obtained with 20, 40, 80, and 160 modes; “direct” indicates the direct integration solution.

Figure 8. Steady-state response of a single beam with a harmonic point force.

Figure 9. Energy densities of the beam at (a) P1, (b) P2, (c) P3, (d) P4, and (e) P5. The vertical dotted lines indicate the time it takes for the wave energy to travel from the center to specified locations.

Figure 10. Energy density distributions at 0 s, 1.2 ms, 2.4 ms, 4.8 ms, 8 ms, and 20 ms (each layer, from top to bottom, is an assembly of results at the specified time).

Figure 11. Coupled beams pinned at both ends.

Figure 12. Steady-state response of the coupled beams.

Figure 13. Energy densities of the coupled beams at (a) P1, (b) P2, (c)P3, (d) P4, and (e) P5 under unloading conditions. The vertical dotted lines indicate the time it takes for the wave energy to travel from the center to specified locations.

Figure 14. Energy densities of the coupled beams at (a) P1, (b) P2, (c)P3, (d) P4, and (e) P5 with a constant input power. The vertical dotted lines indicate the time it takes for the wave energy to travel from the center to specified locations.

Figure 15. Average energy densities of (a) the first beam and (b) the second beam.

Figure 16. Energy densities of the coupled beams at (a) P1, (b) P2, (c) P3, (d) P4, and (e) P5 with a quasi-static input power.

Figure 17. Average energy densities of (a) the first beam and (b) the second beam with a quasi-static input power.

Appendix A

The energy density governing Equation (25) can be solved numerically using FEA. Let the energy density approximate solution of the beam element be of the following form:(A1)$$E(x,t)=N_1(x)\cdot E_1+N_2(x)\cdot E_2={\displaystyle \sum _{i=1}^2{N}_i(x)\cdot E_i}.$$

Choose a collection of liner interpolation functions as follows:(A2)$$N_1(x)=\frac{L-x}{L},N_2(x)=\frac{x}{L}.$$

By applying the Galerkin method of weighted residuals, the residual equation for one beam element can be expressed as:(A3)$${\displaystyle {\int }_0^L\left(\frac{1}{\eta \omega }\frac{{\partial }^{2}E}{\partial t^2}+2\frac{\partial E}{\partial t}+\eta \omega E-\frac{C_g^2}{\eta \omega }{\nabla }^{2}E-{\pi }_{in}\right)\cdot N_j(x)}dx=0,(j=1,2).$$

Substituting Equation (A1) into Equation (A3), the integral items in Equation (A3) can be rearranged as:(A4)$$\begin{array}{c}{\displaystyle {\int }_0^L\frac{1}{\eta \omega }\frac{{\partial }^{2}E}{\partial t^2}\cdot N_j(x)}dx=\frac{1}{\eta \omega }{\displaystyle \sum _{i=1}^2\left({\displaystyle {\int }_0^L{N}_i(x)\cdot N_j(x)}dx\cdot \frac{{\partial }^2{E}_i}{\partial t^2}\right)}\\ {\displaystyle {\int }_0^L\frac{\partial E}{\partial t}\cdot N_j(x)}dx={\displaystyle \sum _{i=1}^2\left({\displaystyle {\int }_0^L{N}_i(x)\cdot N_j(x)}dx\cdot \frac{\partial E_i}{\partial t}\right)}\\ {\displaystyle {\int }_0^L\eta \omega E\cdot N_j(x)}dx=\eta \omega {\displaystyle \sum _{i=1}^2\left({\displaystyle {\int }_0^L{N}_i(x)\cdot N_j(x)}dx\cdot E_i\right)}\\ {\displaystyle {\int }_0^L\frac{C_g^2}{\eta \omega }{\nabla }^{2}E\cdot N_j(x)}dx={\left.\frac{C_g^2}{\eta \omega }\frac{\partial E}{\partial x}\cdot N_j(x)\right|}_0^L-{\displaystyle {\int }_0^L\frac{C_g^2}{\eta \omega }\frac{\partial E}{\partial x}\cdot \frac{\partial N_j(x)}{\partial x}}dx\\ ={\left.-I_j\cdot N_j(x)\right|}_0^L-\frac{C_g^2}{\eta \omega }{\displaystyle \sum _{i=1}^2\left({\displaystyle {\int }_0^L\frac{\partial N_i(x)}{\partial x}\cdot \frac{\partial N_j(x)}{\partial x}}dx\cdot E_i\right)}\end{array}.$$

Substitute Equation (A4) into Equation (A3) and rearrange to yield the FEA-like energy governing equation as follows:(A5)$$M^e\frac{{\partial }^2{E}^e}{\partial t^2}+C^e\frac{\partial E^e}{\partial t}+K^e\cdot E^e=P^e+Q^e,$$ where $M^e$, $C^e$, $K^e$, and ${\pi }^e+Q^e$ correspond to the mass, damping, stiffness, and load matrices, respectively. Their component items are as follows:(A6)$$\begin{array}{c}{M}_{ij}=\frac{1}{\eta \omega }{\displaystyle {\int }_0^L{N}_i(x)\cdot N_j(x)}dx,\\ C_{ij}=2{\displaystyle {\int }_0^L{N}_i(x)\cdot N_j(x)}dx,\\ K_{ij}=\eta \omega {\displaystyle {\int }_0^L{N}_i(x)\cdot N_j(x)}dx+\frac{C_g^2}{\eta \omega }{\displaystyle {\int }_0^L\frac{\partial N_i(x)}{\partial x}\cdot \frac{\partial N_j(x)}{\partial x}}dx,\\ Q_i=-{\left.I_i\cdot N_i(x)\right|}_0^L,\\ P_j={\displaystyle {\int }_0^L{\pi }_{in}\cdot N_j(x)}dx\end{array}.$$

The superscript in Equation (A5) indicates that the equation is for the beam element. After evaluation of the integral terms, the matrices in Equation (A5) can be obtained as:(A7)$$\begin{array}{l}{M}^e=\frac{L}{6\eta \omega }\left[\begin{array}{cc}2& 1\\ 1& 2\end{array}\right],C^e=\frac{L}{3}\left[\begin{array}{cc}2& 1\\ 1& 2\end{array}\right],K^e=\left[\begin{array}{cc}\frac{\eta \omega L}{3}+\frac{C_g^2}{\eta \omega }\cdot \frac{1}{L}& \frac{\eta \omega L}{6}-\frac{C_g^2}{\eta \omega }\cdot \frac{1}{L}\\ \frac{\eta \omega L}{6}-\frac{C_g^2}{\eta \omega }\cdot \frac{1}{L}& \frac{\eta \omega L}{3}+\frac{C_g^2}{\eta \omega }\cdot \frac{1}{L}\end{array}\right],\\ Q^e=\left[\begin{array}{c}{I}_1\\ -I_2\end{array}\right]\end{array}.$$ where $L$ is the element length of the beam.

$Q^e$ in Equation (A7) is the energy intensity that can be cancelled after assembling to system matrices, if no energy reflection or loss occurs at the element boundary. At the system boundary, it is determined by the boundary conditions. Otherwise, a couple matrix $J{C}^e$ is introduced and the energy intensity can be obtained by the energy transmission coefficients [24]:(A8)$$Q_J^e=\left[\begin{array}{cc}1-{\tau }_{11}& -{\tau }_{21}\\ -{\tau }_{12}& 1-{\tau }_{22}\end{array}\right]{\left[\begin{array}{cc}1+{\tau }_{11}& {\tau }_{21}\\ {\tau }_{12}& 1+{\tau }_{22}\end{array}\right]}^{-1}\left\{\begin{array}{c}{C}_{g1}{e}_1\\ C_{g2}{e}_2\end{array}\right\}=J{C}^e\left\{\begin{array}{c}{e}_1\\ e_2\end{array}\right\},$$ where $J{C}^e$ can be expressed as:(A9)$$J{C}^e=\frac{1}{2-{\tau }_{12}-{\tau }_{21}}\left[\begin{array}{cc}{\tau }_{12}& -{\tau }_{21}\\ -{\tau }_{12}& {\tau }_{21}\end{array}\right]\left[\begin{array}{cc}{C}_{g1}& 0\\ 0& C_{g2}\end{array}\right].$$

Assembling the element matrices Equations (A7) and (A9), the transient energy equation can be obtained as:(A10)$$M\frac{{\partial }^{2}E}{\partial t^2}+C\frac{\partial E}{\partial t}+\left(K-JC\right)\cdot E=P_{in},$$ where M, C, and K are constructed from the element matrices Equation (A7). JC is the coupling matrix of the entire system. E is the energy density vector of each node.

The sTEFEA equation can also be derived from Equation (A10) as:(A11)$$\frac{C}{2}\frac{\partial E}{\partial t}+\left(K-JC\right)\cdot E=P_{in}.$$

Appendix B

In this section, the relationship between the energy density and energy intensity is discussed in the quasi-steady state. Taking a two-dimensional plane wave as example, the transverse displacement can be written as:(A12)$$w\left(x,y,t\right)=\left(A_x\cdot e^{-i{K}_x\cdot x}\right)\cdot \left(A_y\cdot e^{-i{K}_y\cdot y}\right)e^{i\omega t},$$ where (A13)$$\begin{array}{l}{K}_x=k_{x1}+k_{x2}=\left(1-\frac{\eta }{4}\right)k_{x1}\\ K_y=k_{y1}+k_{y2}=\left(1-\frac{\eta }{4}\right)k_{y1}\end{array}.$$

The energy density and energy intensity can be represented by the displacement as follows:(A14)$$\begin{array}{c}E=\frac{\rho h}{4}\frac{\partial w}{\partial t}{\left(\frac{\partial w}{\partial t}\right)}^{*}+\\ \frac{D}{4}\left[\frac{{\partial }^{2}w}{\partial x^2}{\left(\frac{{\partial }^{2}w}{\partial x^2}\right)}^{*}+\frac{{\partial }^{2}w}{\partial y^2}{\left(\frac{{\partial }^{2}w}{\partial y^2}\right)}^{*}+2\frac{{\partial }^{2}w}{\partial x^2}{\left(\frac{{\partial }^{2}w}{\partial y^2}\right)}^{*}+2\left(1-\mu \right)\frac{{\partial }^{2}w}{\partial x\partial y}{\left(\frac{{\partial }^{2}w}{\partial x\partial y}\right)}^{*}\right]\\ =\frac{1}{4}\left[\rho h{\omega }^2+D\left({\left|K_x\right|}^4+{\left|K_y\right|}^4+2K_x^2{K}_y^{*2}+2\left(1-\mu \right){\left|K_x\right|}^2{\left|K_y\right|}^2\right)\right]\cdot {\left|A_x\right|}^2{\left|A_y\right|}^2\cdot H\left(x,y\right)\end{array}$$ (A15)$$\begin{array}{c}{I}_x=-\frac{D}{2}\left[\left(\frac{{\partial }^{2}w}{\partial x^2}+\mu \frac{{\partial }^{2}w}{\partial y^2}\right){\left(\frac{{\partial }^{2}w}{\partial x\partial t}\right)}^{*}+\left(1-\mu \right)\frac{{\partial }^{2}w}{\partial x\partial y}{\left(\frac{{\partial }^{2}w}{\partial y\partial t}\right)}^{*}-\frac{\partial }{\partial x}\left(\frac{{\partial }^{2}w}{\partial x^2}+\frac{{\partial }^{2}w}{\partial y^2}\right){\left(\frac{\partial w}{\partial t}\right)}^{*}\right]\\ =\frac{D\omega }{2}\left[K_x^{*}\left(K_x^2+\mu K_y^2\right)+\left(1-\mu \right)K_x{K}_y{K}_y^{*}+K_x\left(K_x^2+K_y^2\right)\right]\cdot {\left|A_x\right|}^2{\left|A_y\right|}^2\cdot H\left(x,y\right)\\ I_y=-\frac{D}{2}\left[\left(\frac{{\partial }^{2}w}{\partial y^2}+\mu \frac{{\partial }^{2}w}{\partial x^2}\right){\left(\frac{{\partial }^{2}w}{\partial y\partial t}\right)}^{*}+\left(1-\mu \right)\frac{{\partial }^{2}w}{\partial x\partial y}{\left(\frac{{\partial }^{2}w}{\partial x\partial t}\right)}^{*}-\frac{\partial }{\partial y}\left(\frac{{\partial }^{2}w}{\partial x^2}+\frac{{\partial }^{2}w}{\partial y^2}\right){\left(\frac{\partial w}{\partial t}\right)}^{*}\right]\\ =\frac{D\omega }{2}\left[K_y^{*}\left(K_y^2+\mu K_x^2\right)+\left(1-\mu \right)K_x{K}_y{K}_x^{*}+K_y\left(K_x^2+K_y^2\right)\right]\cdot {\left|A_x\right|}^2{\left|A_y\right|}^2\cdot H\left(x,y\right)\end{array}$$ where (A16)$$H\left(x,y\right)=e^{-\frac{\eta k_{x1}}{2}\cdot x}\cdot e^{-\frac{\eta k_{y1}}{2}\cdot y}.$$

If assumption (i) is satisfied, the imaginary parts of $K_x$ and $K_y$ can be ignored. Considering that $k^2=k_x^2+k_y^2$, Equations (A14) and (A15) can be rewritten as:(A17)$$\begin{array}{l}{I}_x=D\omega k_x{k}^2\cdot {\left|A_x\right|}^2{\left|A_y\right|}^2\cdot H\left(x,y\right)\\ I_y=D\omega k_y{k}^2\cdot {\left|A_x\right|}^2{\left|A_y\right|}^2\cdot H\left(x,y\right)\\ E=\frac{1}{2}\rho h{\omega }^2\cdot {\left|A_x\right|}^2{\left|A_y\right|}^2\cdot H\left(x,y\right)=\frac{D{k}^4}{2}\cdot {\left|A_x\right|}^2{\left|A_y\right|}^2\cdot H\left(x,y\right)\end{array}.$$

From Equation (A17), the relationship between the energy density and energy intensity can be expressed as:(A18)$$\begin{array}{l}{I}_x=\frac{2\omega k_{x1}}{k^2}E=C_g\cdot E\cdot \cos \theta \\ I_y=\frac{2\omega k_{y1}}{k^2}E=C_g\cdot E\cdot \sin \theta \\ \overrightarrow{I}=\left(I_x,I_y\right)=C_g\cdot E\cdot \left(\cos \theta ,\sin \theta \right)=\overrightarrow{C_g}\cdot E\end{array},$$ where $\theta$ is the angle between the energy intensity and the positive direction of the x-axis. The energy density at each axis is decomposed as follows:(A19)$$\begin{array}{l}{E}_x=E\cdot {\cos }^2\theta ,\\ E_y=E\cdot {\sin }^2\theta ,\\ E=E_x+E_y\end{array}.$$

Then, Equation (A18) can be rewritten as:(A20)$$\begin{array}{l}{I}_x=\frac{C_g}{\cos \theta }\cdot E_x\\ I_y=\frac{C_g}{\sin \theta }\cdot E_y\end{array}.$$

If the reflected waves are considered, the transverse displacement can be written as:(A21)$$w\left(x,y,t\right)=\left(A_x\cdot e^{-i{K}_x\cdot x}+B_x\cdot e^{i{K}_x\cdot x}\right)\cdot \left(A_y\cdot e^{-i{K}_y\cdot y}+B_y\cdot e^{i{K}_y\cdot y}\right)e^{i\omega t}.$$

The energy density and energy intensity can be found in [18]. The energy density and energy intensity should be averaged over one wavelength. If variable substitutions are used as in [17], the energy density can be written as:(A22)$$⟨E⟩=⟨E^{++}⟩+⟨E^{+-}⟩+⟨E^{-+}⟩+⟨E^{--}⟩,$$ where (A23)$$\begin{array}{l}⟨E^{--}⟩=\frac{D{k}^4}{2}\cdot {\left|A_x\right|}^2{\left|A_y\right|}^2\cdot H\left(x,y\right)\\ ⟨E^{-+}⟩=\frac{D{k}^4}{2}\cdot {\left|A_x\right|}^2{\left|B_y\right|}^2\cdot H\left(x,-y\right)\\ ⟨E^{+-}⟩=\frac{D{k}^4}{2}\cdot {\left|B_x\right|}^2{\left|A_y\right|}^2\cdot H\left(-x,y\right)\\ ⟨E^{++}⟩=\frac{D{k}^4}{2}\cdot {\left|B_x\right|}^2{\left|B_y\right|}^2\cdot H\left(-x,-y\right)\end{array}.$$

The energy density of each axis is decomposed as follows:(A24)$$\begin{array}{l}⟨E_x^{+}⟩=\left(⟨E^{-+}⟩+⟨E^{--}⟩\right)\cdot {\cos }^2\theta \\ ⟨E_x^{-}⟩=\left(⟨E^{++}⟩+⟨E^{+-}⟩\right)\cdot {\cos }^2\theta \\ ⟨E_y^{+}⟩=\left(⟨E^{+-}⟩+⟨E^{--}⟩\right)\cdot {\sin }^2\theta \\ ⟨E_y^{-}⟩=\left(⟨E^{++}⟩+⟨E^{-+}⟩\right)\cdot {\sin }^2\theta \end{array},$$ and (A25)$$\begin{array}{l}⟨E_x⟩=⟨E_x^{+}⟩+⟨E_x^{-}⟩=E\cdot {\cos }^2\theta \\ ⟨E_y⟩=⟨E_y^{+}⟩+⟨E_y^{-}⟩=E\cdot {\sin }^2\theta \end{array}$$

The total energy density can be rewritten as:(A26)$$⟨E⟩=⟨E_x^{+}⟩+⟨E_x^{-}⟩+⟨E_y^{+}⟩+⟨E_y^{-}⟩.$$

The energy intensity components in each direction are as follows:(A27)$$\begin{array}{l}⟨I_x^{+}⟩=D\omega k_x{k}^2\cdot \left({\left|A_x\right|}^2{\left|A_y\right|}^2\cdot H\left(x,y\right)+{\left|A_x\right|}^2{\left|B_y\right|}^2\cdot H\left(x,-y\right)\right)\\ ⟨I_x^{-}⟩=D\omega k_x{k}^2\cdot \left({\left|B_x\right|}^2{\left|A_y\right|}^2\cdot H\left(-x,y\right)+{\left|B_x\right|}^2{\left|B_y\right|}^2\cdot H\left(-x,-y\right)\right)\\ ⟨I_y^{+}⟩=D\omega k_y{k}^2\cdot \left({\left|A_x\right|}^2{\left|A_y\right|}^2\cdot H\left(x,y\right)+{\left|B_x\right|}^2{\left|A_y\right|}^2\cdot H\left(-x,y\right)\right)\\ ⟨I_y^{-}⟩=D\omega k_y{k}^2\cdot \left({\left|A_x\right|}^2{\left|B_y\right|}^2\cdot H\left(x,-y\right)+{\left|B_x\right|}^2{\left|B_y\right|}^2\cdot H\left(-x,-y\right)\right)\end{array}.$$

Comparing Equations (A24) and (A27), a similar relationship to Equation (A18) can be expressed as:(A28)$$\begin{array}{l}⟨I_x^{\pm }⟩=\pm C_g\cdot \cos \theta \cdot \left(⟨E^{\pm +}⟩+⟨E^{\pm -}⟩\right)=\pm \frac{C_g}{\cos \theta }\cdot ⟨E_x^{\pm }⟩\\ ⟨I_y^{\pm }⟩=\pm C_g\cdot \sin \theta \cdot \left(⟨E^{+\pm }⟩+⟨E^{-\pm }⟩\right)=\pm \frac{C_g}{\sin \theta }\cdot ⟨E_y^{\pm }⟩\end{array}.$$

The energy intensity vector can be expressed as:(A29)$$\overrightarrow{I}=\left(I_x^{-}+I_x^{+},I_y^{-}+I_y^{+}\right)=\left(I_x,I_y\right).$$

If only the steady state is considered, the relationship between the energy density and energy intensity can be expressed as:(A30)$$\left\{\begin{array}{l}⟨I_i^{+}⟩=-\frac{C_g^2}{\eta \omega }\cdot \frac{\partial ⟨E_i^{+}⟩}{\partial x_i}\\ ⟨I_i^{-}⟩=\frac{C_g^2}{\eta \omega }\cdot \frac{\partial ⟨E_i^{-}⟩}{\partial x_i}\\ ⟨I_i⟩=⟨I_i^{+}⟩+⟨I_i^{-}⟩\end{array}\right.,\left(i=x,y\right).$$

References

1. Lyon, R.H.; DeJong, R.G. Theory and Application of Statistical Energy Analysis; 2nd ed. Newnes: Boston, MA, USA, 1995; [DOI: https://dx.doi.org/10.1016/C2009-0-26747-X]

2. Desmet, W. Mid-frequency vibro-acoustic modelling: Challenges and potential solutions. Proceedings of the ISMA; Lueven, Belgium, 16–18 September 2002; 12.

3. Roque, C.M.C.; Ferreira, A.J.M.; Reddy, J.N. Analysis of Timoshenko nanobeams with a nonlocal formulation and meshless method. Int. J. Eng. Sci.; 2011; 49, pp. 976-984. [DOI: https://dx.doi.org/10.1016/j.ijengsci.2011.05.010]

4. Banerjee, J.R. Free vibration analysis of a twisted beam using the dynamic stiffness method. Int. J. Solids Struct.; 2001; 38, pp. 6703-6722. [DOI: https://dx.doi.org/10.1016/S0020-7683(01)00119-6]

5. Banerjee, J.R.; Ananthapuvirajah, A. An exact dynamic stiffness matrix for a beam incorporating Rayleigh–Love and Timoshenko theories. Int. J. Mech. Sci.; 2019; 150, pp. 337-347. [DOI: https://dx.doi.org/10.1016/j.ijmecsci.2018.10.012]

6. Wang, T.; Altabey, W.A.; Noori, M.; Ghiasi, R. A Deep Learning Based Approach for Response Prediction of Beam-like Structures. Struct. Durab. Health Monit.; 2020; 14, 315. [DOI: https://dx.doi.org/10.32604/sdhm.2020.010564]

7. Rajasekaran, S.; Khaniki, H.B.; Ghayesh, M.H. On the mechanics of shear deformable micro beams under thermo-mechanical loads using finite element analysis and deep learning neural network. Mech. Based Des. Struct. Mach.; 2022; 51, pp. 6612-6656. [DOI: https://dx.doi.org/10.1080/15397734.2022.2047721]

8. Liu, Z.; Niu, J.; Gao, X. An improved approach for analysis of coupled structures in Energy Finite Element Analysis using the coupling loss factor. Comput. Struct.; 2018; 210, pp. 69-86. [DOI: https://dx.doi.org/10.1016/j.compstruc.2018.08.017]

9. Manning, J.E. Predicting mechanical shock transmission. Shock Vib. Bull.; 1968; 37, pp. 65-70.

10. Pinnington, R.J.; Lednik, D. Transient statistical energy analysis of an impulsively excited two oscillator system. J. Sound Vib.; 1996; 189, pp. 249-264. [DOI: https://dx.doi.org/10.1006/jsvi.1996.0018]

11. Pinnington, R.J.; Lednik, D. Transient energy flow between two coupled beams. J. Sound Vib.; 1996; 189, pp. 265-287. [DOI: https://dx.doi.org/10.1006/jsvi.1996.0019]

12. Langley, R.S.; Hawes, D.H.; Butlin, T.; Ishii, Y. A derivation of the Transient Statistical Energy Analysis (TSEA) equations with benchmark applications to plate systems. J. Sound Vib.; 2019; 445, pp. 88-102. [DOI: https://dx.doi.org/10.1016/j.jsv.2019.01.002]

13. Hawes, D.H.; Langley, R.S.; Butlin, T.; Ishii, Y. A hybrid Finite Element-Statistical Energy Analysis method for impulsive and transient loading. J. Sound Vib.; 2019; 459, 114849. [DOI: https://dx.doi.org/10.1016/j.jsv.2019.114849]

14. Nefske, D.J.; Sung, S.H. Power Flow Finite Element Analysis of Dynamic Systems: Basic Theory and Application to Beams. J. Vib. Acoust.; 1989; 111, pp. 94-100. [DOI: https://dx.doi.org/10.1115/1.3269830]

15. Wohlever, J.; Bernhard, R. Energy distributions in rods and beams. Proceedings of the 12th Aeroacoustic Conference; San Antonio, TX, USA, 10–12 April 1989; 1122.

16. Bouthier, O.; Bernhard, R. Simple models of energy flow in vibrating membranes. J. Sound Vib.; 1995; 182, pp. 129-147. [DOI: https://dx.doi.org/10.1006/jsvi.1995.0186]

17. Bouthier, O.M.; Bernhard, R.J. Models of space-averaged energetics of plates. AIAA J.; 1992; 30, pp. 616-623. [DOI: https://dx.doi.org/10.2514/3.10964]

18. Bouthier, O.; Bernhard, R. Simple models of the energetics of transversely vibrating plates. J. Sound Vib.; 1995; 182, pp. 149-164. [DOI: https://dx.doi.org/10.1006/jsvi.1995.0187]

19. Park, D.-H.; Hong, S.-Y.; Kil, H.-G.; Jeon, J.-J. Power flow models and analysis of in-plane waves in finite coupled thin plates. J. Sound Vib.; 2001; 244, pp. 651-668. [DOI: https://dx.doi.org/10.1006/jsvi.2000.3517]

20. Park, Y.-H.; Hong, S.-Y. Vibrational Energy Flow Analysis of Corrected Flexural Waves in Timoshenko Beam—Part I: Theory of an Energetic Model. Shock Vib.; 2006; 13, pp. 137-165. [DOI: https://dx.doi.org/10.1155/2006/308715]

21. Park, Y.-H.; Hong, S.-Y. Vibrational power flow models for transversely vibrating finite Mindlin plate. J. Sound Vib.; 2008; 317, pp. 800-840. [DOI: https://dx.doi.org/10.1016/j.jsv.2008.03.049]

22. Park, Y.-H. Energy flow finite element analysis of general Mindlin plate structures coupled at arbitrary angles. Int. J. Nav. Archit. Ocean Eng.; 2019; 11, pp. 435-447. [DOI: https://dx.doi.org/10.1016/j.ijnaoe.2018.08.001]

23. Han, J.-B.; Hong, S.-Y.; Song, J.-H.; Kwon, H.-W. Vibrational energy flow models for the Rayleigh–Love and Rayleigh–Bishop rods. J. Sound Vib.; 2014; 333, pp. 520-540. [DOI: https://dx.doi.org/10.1016/j.jsv.2013.08.027]

24. Yan, X. Energy Finite Element Analysis Developments for High Frequency Vibration Analysis of Composite Structures. Ph.D. Thesis; The University of Michigan: Ann Arbor, MI, USA, 2008.

25. Lin, Z.; Chen, X.; Zhang, B. Application of the energy finite element analysis to vibration of beams with stepped thickness and variable cross-section. J. Vibroeng.; 2018; 20, pp. 2237-2252. [DOI: https://dx.doi.org/10.21595/jve.2018.19625]

26. Zhang, W.; Vlahopoulos, N.; Wu, K. An energy finite element formulation for high-frequency vibration analysis of externally fluid-loaded cylindrical shells with periodic circumferential stiffeners subjected to axi-symmetric excitation. J. Sound Vib.; 2005; 282, pp. 679-700. [DOI: https://dx.doi.org/10.1016/j.jsv.2004.03.063]

27. Zhang, W.; Wang, A.; Vlahopoulos, N.; Wu, K. A vibration analysis of stiffened plates under heavy fluid loading by an energy finite element analysis formulation. Finite Elem. Anal. Des.; 2005; 41, pp. 1056-1078. [DOI: https://dx.doi.org/10.1016/j.finel.2004.10.012]

28. Xie, M.; Li, L.; Shang, X.; Zhao, J. Initial Investigation of Energy Finite Element Validation on High-Frequency Flexural Vibration of Stiffened Thin Orthotropic Plates. Shock Vib.; 2018; 2018, 10. [DOI: https://dx.doi.org/10.1155/2018/1461086]

29. Nokhbatolfoghahai, A.; Navazi, H.; Haddadpour, H. High-frequency random vibrations of a stiffened plate with a cutout using energy finite element and experimental methods. Proc. Inst. Mech. Eng. Part C J. Mech. Eng. Sci.; 2020; 234, pp. 3297-3317. [DOI: https://dx.doi.org/10.1177/0954406220914328]

30. Teng, X.; Han, Y.; Jiang, X.; Chen, X.; Zhou, M. Energy Flow Analysis Model of High-Frequency Vibration Response for Plates with Free Layer Damping Treatment. Mathematics; 2023; 11, 1379. [DOI: https://dx.doi.org/10.3390/math11061379]

31. Wang, D.; Xie, M.; Li, Y. High-Frequency Dynamic Analysis of Plates in Thermal Environments Based on Energy Finite Element Method. Shock Vib.; 2015; 2015, 14. [DOI: https://dx.doi.org/10.1155/2015/157208]

32. Chen, Z.; Yang, Z.; Guo, N.; Zhang, G. An energy finite element method for high frequency vibration analysis of beams with axial force. Appl. Math. Model.; 2018; 61, pp. 521-539. [DOI: https://dx.doi.org/10.1016/j.apm.2018.04.016]

33. Xie, M.X.; Chen, H.L.; Wu, J.H.; Sun, F.G. Application of Energy Finite Element Method to High-frequency Structural-acoustic Coupling of an Aircraft Cabin with Truncated Conical Shape. Comput. Model. Eng. Sci.; 2010; 61, pp. 1-22. [DOI: https://dx.doi.org/10.3970/cmes.2010.061.001]

34. Vlahopoulos, N.; Schiller, N.; Lee, S. Energy finite element analysis developments for vibration analysis of composite aircraft structures. SAE Int. J. Aerosp.; 2011; 4, pp. 593-601. [DOI: https://dx.doi.org/10.4271/2011-01-1734]

35. Shen, Y.; Xu, Y.; Sheng, X.; Wan, S.; Yin, X. Vibration Energy Transfer Characteristics of Panels with Multiple Coupling Forms in Satellites. Int. J. Aeronaut. Space Sci.; 2023; 24, pp. 1231-1243. [DOI: https://dx.doi.org/10.1007/s42405-023-00576-6]

36. Zhang, W.; Wang, A.; Vlahopoulos, N. An alternative energy finite element formulation based on incoherent orthogonal waves and its validation for marine structures. Finite Elem. Anal. Des.; 2002; 38, pp. 1095-1113. [DOI: https://dx.doi.org/10.1016/S0168-874X(02)00053-7]

37. Homg, S.B.; Vlahopoulos, N. Application of a hybrid finite-element formulation for structure-borne vibration in an automotive vehicle structure. J. Acoust. Soc. Am.; 2006; 119, pp. 3297-3298. [DOI: https://dx.doi.org/10.1121/1.4786242]

38. Vlahopoulos, N.; Wu, K. Energy finite element analysis of naval vehicles. J. Acoust. Soc. Am.; 2010; 127, 1850. [DOI: https://dx.doi.org/10.1121/1.3384371]

39. Zheng, X.; Dai, W.; Qiu, Y.; Hao, Z. Prediction and energy contribution analysis of interior noise in a high-speed train based on modified energy finite element analysis. Mech. Syst. Signal Process.; 2019; 126, pp. 439-457. [DOI: https://dx.doi.org/10.1016/j.ymssp.2019.02.042]

40. Ichchou, M.N.; Le Bot, A.; Jezequel, L. A transient local energy approach as an alternative to transient sea: Wave and telegraph equations. J. Sound Vib.; 2001; 246, pp. 829-840. [DOI: https://dx.doi.org/10.1006/jsvi.2001.3647]

41. Sui, F.S.; Ichchou, M.N.; Jezequel, L. Prediction of vibroacoustics energy using a discretized transient local energy approach and comparison with TSEA. J. Sound Vib.; 2002; 251, pp. 163-180. [DOI: https://dx.doi.org/10.1006/jsvi.2001.3786]

42. Sui, F.; Ichchou, M.N. Prediction of Time-Varying Vibroacoustic Energy Using a New Energy Approach. J. Vib. Acoust.; 2004; 126, pp. 184-189. [DOI: https://dx.doi.org/10.1115/1.1687399]

43. Chen, Q.; Fei, Q.; Li, Y.; Wu, S.; Yang, X. Prediction of the transient energy response for complex vibro-acoustic systems. J. Mech. Sci. Technol.; 2019; 33, pp. 495-504. [DOI: https://dx.doi.org/10.1007/s12206-019-0102-y]

44. Cremer, L.; Heckl, M.; Petersson, B.A.T. Structure-Borne Sound: Structural Vibrations and Sound Radiation at Audio Frequencies; 3rd ed. Springer: Berlin/Heidelberg, Germany, 2005.

45. Cho, P.E.-H. Energy Flow Analysis of Coupled Structures. Ph.D. Thesis; Purdue University: West Lafayette, IN, USA, 1993.

46. Langley, R.S.; Heron, K.H. Elastic wave transmission through plate/beam junctions. J. Sound Vib.; 1990; 143, pp. 241-253. [DOI: https://dx.doi.org/10.1016/0022-460X(90)90953-W]

47. Lai, M.; Soom, A. Statistical energy analysis for the time-integrated transient response of vibrating systems. J. Vib. Acoust.; 1990; 112, pp. 206-213. [DOI: https://dx.doi.org/10.1115/1.2930114]

48. Kong, X.; Chen, H.; Zhu, D.; Zhang, W. Study on the validity region of Energy Finite Element Analysis. J. Sound Vib.; 2014; 333, pp. 2601-2616. [DOI: https://dx.doi.org/10.1016/j.jsv.2013.11.047]

49. Chen, Q.; Fei, Q.; Wu, S.; Li, Y.; Yang, X. An efficient transient analysis method for time-varying structures based on statistical energy analysis. Mech. Res. Commun.; 2018; 91, pp. 93-99. [DOI: https://dx.doi.org/10.1016/j.mechrescom.2018.06.001]