Energy flow boundary element method for vibration analysis of one and two dimension structures

In this paper, Energy Flow Boundary Element Method (EFBEM) was developed to predict the vibration behavior of oneand two-dimensional structures in the medium-to-high frequency ranges. Free Space Green functions used in the method were obtained from EFA energy equations. Direct and indirect EFBEMs were formulated for both oneand two-dimensional cases, and numerically applied to predict the energy density and intensity distributions of simple Euler-Bernoulli beams, single rectangular thin plates, and L-shaped thin plates vibrating in the medium-to-high frequency ranges. The results from these methods were compared with the EFA solutions to verify the EFBEM.


List of symbols
∂: partial derivative t: time ∇: gradient q: intensity vector q: intensity π: power < >: time-average η: structural damping loss factor ω: circular frequency e: total energy density c g : group velocity δ(x): Dirac delta function G: free space Green's function F : fundamental solution φ: fictitious sources Ω: domain Γ: boundary n: normal vector E: Young's modulus ρ: mass density f : frequency

Introduction
Many noise and vibration engineers and investigators have shown much interest in the dynamic complex structures vibrating in the medium-to-high frequency ranges.The traditional finite element method (FEM) and boundary element method (BEM) have been widely used for the low frequency analysis.At high frequencies, these methods require a higher order shape function or more fine elements to obtain more accurate results.Thus, the traditional FEM and BEM have become costly and ineffective for vibration analysis of complex structures in high frequency ranges.The statistical energy analysis (SEA), which has been developed for high frequency vibration analysis, gives single averaged value for the energy density of a subsystem, and thus it cannot provide detailed information such as variation of energy density and energy flow paths in a subsystem, and SEA requires that the size of subsystem in a structure be large when the frequency of interest is decreased.
To improve the weaknesses of SEA in vibration analysis in the high frequency ranges and to overcome the frequency-dependent limitations of the traditional FEM and BEM, researchers have investigated alternative methods.One of the promising alternative method is Energy Flow Analysis (EFA), which was introduced by Belov et al. [1] in 1977.EFA method is based on an energy governing equation analogous to the heat conductivity equation, in which the main variable is the energy density.Using this method at high frequencies, the spatial variation of the timeand locally space-averaged vibration energy density and energy transmission paths in a structure can be effectively predicted.Moreover, unlike the traditional FEM, BEM and SEA, the EFA method has been considered appropriate for the vibration analysis of complex structures in mid-frequency ranges.Nefske and Sung [2] implemented the finite element formulation of the energy governing equation of EFA, and used it to predict the vibration responses of an uncoupled beam and two coupled beams excited by a harmonic force.Wohlever and Bernhard [3,4] derived the energy governing equations for the propagation of the vibration waves in rods and Euler-Bernoulli beams, and Bouthier and Bernhard [4][5][6][7] derived the energy governing equations for the flexural waves in membranes and thin plates and for the acoustic waves in enclosures.Cho [8] studied the joint relationships of coupled structures by using the concept of power transmission and reflection coefficients.Park et al. [9] developed the energy flow models of the in-plane waves in isotropic thin plates and the flexural waves in orthotropic thin plates, and Seo et al. [10] researched the energy flow analysis of beam-plate coupled structures.Similar to EFA, the simplified energy method (SEM) was developed by Lase, Le Bot et al. [11,12]; it is a simplification of the general energy method (GEM).
This work aims to implement boundary element formulations of the energy governing equation of EFA for oneand two-dimensional waves [13,14].Direct and indirect methods are investigated and numerically applied to the prediction of the vibration energy density and intensity distributions in simple beams, rectangular thin plates, and L-type thin plates.

Energy governing equation
In steady state (∂ / ∂ t = 0), the vibrational power injected into a structural element is equal to the sum of the power flowing out through its boundary and the power dissipated in the element.The steady state energy balance equation of an elastic structural element has the form of where q is the intensity vector, which means the power transmitted through the boundary of the element, and π diss is the dissipated power due to the structural damping and π in is the power injected due to external loads.The timeaveraged dissipated power in an elastic medium with small structural damping is proportional to the time-averaged total energy density, and this concept can be expressed as where the brackets < > mean the time-average over a period, η is the structural damping loss factor, ω is the circular frequency, and e is the total energy density.In Eq. ( 2), it is assumed that the kinetic energy and the potential energy of the medium are approximately the same.
Based on the work of Wohlever and Bernhard, the time-averaged intensity of longitudinal and torsional waves in rods can be expressed as the derivative of the total energy density as follows: where c g is the group velocity of the corresponding wave.The intensity and energy density of flexural waves in Euler-Bernoulli beams have the same energy transmission relation as Eq. ( 3).In this case, < q > and < e > in Eq. ( 3) represent the time-and locally space-averaged intensity and energy density of far-field components, respectively.The combination of Eqs (1), ( 2) and (3) yields the energy governing equation of EFA for one-dimensional waves in the form of When two-dimensional in-plane waves in thin plates and flexural waves in membranes are considered, the energy transmission relation is more generally expressed as and thus, the energy governing equation is written by where the brackets < >mean the time-and space-average over a period and a half wavelength.When flexural waves propagate in thin plates, the time-and locally space-averaged intensity and energy density of far-field components satisfy Eqs ( 5) and (6).

Formulation of the direct method
In this section, the direct boundary element method is applied to solve the energy governing equations for the vibrational waves in one-dimensional structures.First, Eq. ( 4), the energy governing equation for one-dimensional structures, can be rewritten implicitly as where the brackets < > are neglected for brevity.If Eq. ( 7) is multiplied by a two-variable function G(x; ξ) and is integrated over the range 0 x L, a weak variational form of Eq. ( 4) can be obtained by where L is the length of the structure, and G(x; ξ) is specified to be the solution of the equation where δ(x − ξ) is the Dirac delta function, which is mathematically equivalent to the effect of a unit concentrated source applied at a point ξ.When the function G(x; ξ) is the fundamental solution (free space Green's function) of Eq. ( 9), we can obtain the following expression: If Eq. ( 8) is integrated by parts twice, we obtain and when Eqs (3) and ( 9) are substituted into Eq.( 11), the energy density at a point ξ in the one dimensional domain can be obtained by where F (x; ξ) is defined as and ξ is the field point coordinate in the domain and x is the source point on the boundary.By using Eq. ( 3) and taking the first derivative of e(ξ) with respect to the field variable ξ from Eq. ( 12), the intensity is written as follows:

Formulation of the indirect method
A key concept of the indirect boundary element method is the embedding of the real system in the infinite field to produce a fictitious system, assuming that a fictitious source exists on the original boundary of the real system.When this assumption is applied to solve Eq. ( 4), the energy governing equation for the vibrational problems of one-dimensional systems in medium-to-high frequency ranges, and the vibrational energy density in the domain can be expressed as follows: and from Eqs (3) and (15), we can obtain the intensity such as where x is the field point coordinate in the domain contrary to the direct method.G(x; ξ) and F (x; ξ) are the fundamental solution and its first derivative, as in Eqs ( 10) and ( 13), respectively.Also, φ and π in are the symbols for sources, and the symbols are distinct in that π in is used for known, specified source strengths, whereas φ is defined solely for fictitious sources, φ(0) and φ(L), applied at the boundaries of a fictitious system.If the energy density is given for the boundary condition, Eq. ( 15) is used; if the intensity given for the boundary condition, Eq. ( 16) is used.And so φ(0) and φ(L) are calculated first on the boundary of the system.Next, we can obtain the distributions of the energy density and intensity in the domain when these sources are substituted into Eqs (15) and ( 16) and these equations are solved.

Formulation of the direct method
In this section, to analyze the vibrational phenomena for two dimensional structures such as a thin plate and a membrane in medium-to-high frequency ranges, the direct boundary element method is applied to Eq. ( 6) in a way similar to one dimension formulation presented at Section 3.1.
Using the weighted residual technique, Eq. ( 6) can be expressed implicitly as If Eq. ( 17) is multiplied by a weight function G( x, ξ) and is integrated over the interested domain Ω in two dimensional structures, we can obtain the following equation where G( x, ξ) denotes the solution of the energy governing equation with a unit input power as follows: and is expressed by where G( x, ξ) is the fundamental solution for an isotropic two-dimensional region.Here, K 0 is the zeroth order modified Bessel function of the second kind and the distance from the field point ξ to the source point With the application of Green's second identities providing a relationship between the surface integral and the line integral over the line bounding the original domain, Eq. ( 21) can be rewritten as where F ( x, ξ) is defined as and n( x) denotes the components of a unit vector at the source point x defining the outward normal direction to the boundary Γ. q n is the outward normal intensity to the boundary and is expressed by Additionally, the coefficient c( ξ) is α/2π in a two-dimensional case, where α means the internal angle (in radian) at a field point ξ.22) indicates the vibrational energy density at a point ξ, the intensity of the ξ i -direction in the domain is obtained by taking the first derivative of Eq. ( 22) with respect to the field variable ξ i as follows: where ξ i (for i = 1, 2) is the components of a rectangular Cartesian coordinate system.

Formulation of the indirect method
To keep the basic concept of the indirect method, we assume that the real system is embedded in an infinite two-dimensional structure, and a fictitious source is distributed on the boundary line.With this assumption, the energy density and the outward normal intensity to the boundary can be expressed by and respectively.x is a field point in Ω, and z is the position vector of the input power π in .Thus, the fictitious source φ( ξ) means initially the unknown intensity per unit length of Γ at a source point ξ on Γ. G( x, ξ) and F ( x, → ξ ) are identical to Eqs (20) and (24).

Numerical application for the one dimensional problem
In this section, the energy density and the intensity distributions for the vibrational waves of a simple Euler-Bernoulli beam structure are predicted by using the one dimensional energy flow boundary element formulations in the form of the direct method and the indirect method, Eqs ( 12), ( 14), ( 15) and ( 16).These equations are used with no discretization because the boundaries of the one-dimensional field are merely two points, and therefore, only two boundary elements are involved.If a simply-supported beam structure is excited by a transverse harmonic point force at a single frequency as shown in Fig. 1(a), zero energy outflow (intensity) boundary conditions at the beam edges can be assumed like Fig. 1(b).The Young's modulus of the beam is E = 7.1 × 10 10 N/m 2 and the mass density is ρ = 2700 kg/m 3 .The length of the beam is L = 12 m and the exciting force is located at x = 4 m from the left end of the beam.
In the first example, the exciting frequency is f = 1 kHz and the damping loss factor of the beam is η = 0.02.and the damping loss factor is η = 0.02, and from these results, we can confirm that EFBEM can be successfully implemented for one-dimensional problems.

Numerical application for the two dimensional problem
When the energy density and intensity of two-dimensional structures in medium-to-high frequency ranges are analyzed, Eqs ( 22) and (25) for the direct EFBEM are used, and Eqs (26) and ( 27) for the indirect EFBEM are applied.If we were able to integrate these equations in closed form and solve them for the fictitious source φ or the initially known, boundary values of e and q, then the solution would be exact.But this operation is virtually impossible for practical problems, and therefore, the domain and boundary have to be discretized, trading computing time for solution accuracy.
This discretization scheme utilizes constant boundary elements, characterized by their midpoints, and uniform distribution of variables (e.g.e, q and φ) over the elements.Discrete approximation of Eq. ( 22) for e( x i 0 ), the energy density on the ith boundary element, in the direct EFBEM is written as where ξ i 0 is the position vector of the midpoint of the ith boundary element, ∆Γ is the length of the jth boundary element, z k is the position vector of the midpoint of the kth internal cell, and ∆Ω is the area of the kth internal cell.Thus, Γ and Ω are approximated by N boundary segments and M internal cells, respectively.The intensity obtained from Eq. ( 25) can be expressed in the same way as Eq.(28).Equations ( 26) and ( 27) of the indirect EFBEM are rewritten as Eqs ( 29) and (30) by discretization: where x i 0 is the position vector of the midpoint of the ith boundary element.If linear elements are used in the discretizaton, discrete approximations similar to Eqs (28), ( 29) and (30) can be easily derived except that e, q n and φ cannot be taken out of the integral sign and c( ξ) is not 1/2.From discrete equations applied to the boundary conditions, the unknown boundary values are determined, and the energy density and intensity in the domain are finally obtained.As the first example, the rectangular thin plate of Fig. 4(a) applied with a transverse harmonic force is considered, and simply-supported edges can be assumed by zero energy outflow (intensity) boundary conditions as in Fig. 4(b).Then, the dimensions and thickness of the plate are L x = L y = 1 m and h = 1 mm, respectively, and the material properties of the plate are those of aluminum.The force is located at x = 0.5 m and y = 0.5 m in the plate and its magnitude is 1 N.
When the frequency is f = 1 kHz and the damping loss factor is η = 0.2, the results of Fig. 5 show the energy density distributions calculated from various methods.Figure 5(a) shows the classical solution [8,9] of the energy density obtained from the motion equation of the classical plate theory, and Fig. 5(b) is the EFA solution [8,9] calculated from Eq. ( 6), the energy governing equation for the transverse vibration of the plate by the double Fourier series.To assure the convergence of the series, the lowest 22,500 modes are calculated and summed.The energy densities predicted by the direct EFBEM and indirect EFBEM are shown in Fig. 5(c) and 5(d), respectively.Near the driving point, generally, EFBEM result has a sharp peak because the EFBEM solution is singular, like 1/ √ r, where r is the source-field distance.Accordingly, the neighbor of the driving point is excluded from field points.Figure 6 shows the comparison of the energy density distributions in a rectangular plate predicted by classical solution, SEA solution, EFA analytic solution, direct EFBEM solution and indirect EFBEM solution along the representative line of Fig. 4(b).The results of both EFBEM are similar to those of other methods.The error rates between EFA analytic solution and both EFBEM solutions in Fig. 6 are shown in Fig. 7.The error rates of both methods are similar in the direct field near the actual source, but the error of the indirect method is more than that of direct method in the diffracted field or scattered field near the secondary source on the boundary.
When the frequency is f = 10 kHz and the damping loss factor is η = 0.05, The error rates of between energy densities calculated by the EFA and direct EFBEM are compared as the number of boundary elements is changed, as shown in Fig. 8.
Like the traditional BEM, when there are more elements, the EFBEM solutions become more accurate, but EFBEM solution calculated with a small number of boundary elements is sufficiently accurate at high frequency.
For a EFBEM model composed of 64 boundary elements, the error rates of energy density at various frequencies for a constant damping loss factor (η = 0.02) are compared in Fig. 9, and the error rates of energy density at various damping loss factors for a constant frequency (f = 5 kHz) are compared in Fig. 10.The pattern of error rate between the EFA solution and direct EFBEM solution varies with the increase of frequency as shown in Fig. 9(a), but the difference between both solutions decreases as the frequency increases, as shown in Fig. 9    EFBEM offers a reliable result at high frequency.The results of Fig. 10(a) and 10(b) show that very light damping loss factor causes a large error; especially, when damping loss factor is 0.001, the error is about 2 × 10 −6 .The error of EFBEM is mainly generated near the source and the boundary, and the error near the source is more than that near the boundary as the frequency and damping loss factor increase.The above results verify the accuracy of EFBEM and confirms that EFBEM can be efficiently applied for vibration problems of two-dimensional structures.
In the second example, energy flow boundary element analysis for an L-shaped thin plate, which cannot be applied to the classical solution, is performed.The L-shaped plate in Fig. 11 is excited by a transverse harmonic point force at location x = 0.8 m and y = 0.25 m with zero energy flow boundary conditions.The plate is 1 m × 1 m × 1 mm in dimensions and has the material properties of aluminum.For the frequency of f = 1 kHz and the damping of η = 0.02, the energy density and the intensity evaluated by the direct EFBEM are shown as in Fig. 12.The direct EFBEM analysis results for the case in which damping loss factor is increased to η = 0.2 at the same exciting frequency are illustrated in Fig. 13.The increase of the damping causes large spatial variation of the energy density, decreases the whole value of the energy density, and also rapidly diminishes the magnitude of the intensity, as shown in Fig. 13

Conclusions
In this work, energy flow boundary element method has been newly developed as a tool for the medium-to-high frequency vibration analysis of one-and two-dimensional structures.To derive the energy integral formulations, free space Green functions obtained from the EFA energy equation is used and the fundamental approaches for the   Further studies are recommended on the development of EFBEM for three-dimensional vibration problems and on the investigation of the acoustic energy density and intensity field using EFBEM.

Fig. 1 .
Fig. 1.Simply-supported beam.(a) Geometry of beam with transverse harmonic load.(b) Energy flow model with dimensions and zero energy flow boundary conditions.

Figure 2 (
a) is the comparison of the results which are the EFA analytic solution of the energy density calculated from the energy governing equation and the numerical solutions of the energy density obtained from the direct method

Fig. 2 .Fig. 3 .
Fig. 2. Comparisons of energy density distributions in beam when η = 0.02.The reference energy density is 1 × 10 −12 J/m 2 ; (a) f = 1 kHz, (b) f = 10 kHz: -, EFA; ---, SEA; , direct EFBEM; , indirect EFBEM.and the indirect method of the energy flow boundary element method.From Fig. 2(a), we know that the results of EFBEM for one-dimensional problems agree fairly well with the EFA analytic solution.The energy density distribution evaluated by the statistical energy analysis (SEA) is compared with that from EFBEM.EFA expresses the spatial variation of the energy density better than SEA.Moreover, EFBEM gives highly accurate values of intensity as shown in Fig. 3(a).Figures 2(b) and 3(b) are the EFBEM results when the frequency is f = 10 kHz

Fig. 4 .
Fig. 4. Simply-supported rectangular plate.(a) Geometry of plate with external transverse harmonic load.(b) Energy flow model with dimensions and zero energy flow boundary conditions.

Fig. 11 .
Fig. 11.Simply-supported L-shaped plate.(a) Geometry of plate with external transverse harmonic load.(b) Energy flow model with dimensions and zero energy flow boundary conditions. (a).
direct and indirect method are performed.The developed energy integral equations are discretized and applied to predict the energy density and intensity distributions of a simple beam and a thin plate.From numerical examples of the beam and plate, it has been confirmed that the approximate energy density and intensity field obtained by energy flow boundary element analysis well represent the global variation of the response.