Vibrational energy flow analysis of corrected flexural waves in Timoshenko beam – Part I : Theory of an energetic model

In this paper, an energy flow model is developed to analyze transverse vibration including the effects of rotatory inertia as well as shear distortion, which are very important in the Timoshenko beam transversely vibrating in the medium-to-high frequency ranges. The energy governing equations for this energy flow model are newly derived by using classical displacement solutions of the flexural motion for the Timoshenko beam, in detail. The derived energy governing equations are in the general form incorporating not only the Euler-Bernoulli beam theory used for the conventional energy flow model but also the Rayleigh, shear, and Timoshenko beam theories. Finally, to verify the validity and accuracy of the derived model, numerical analyses for simple finite Timoshenko beams were performed. The results obtained by the derived energy flow model for simple finite Timoshenko beams are compared with those of the classical solutions for the Timoshenko beam, the energy flow solution, and the classical solution for the Euler-Bernoulli beam with various excitation frequencies and damping loss factors of the beam. In addition, the vibrational energy flow analyses of coupled Timoshenko beams are described in the other companion paper.


Introduction
Recently, the noise and vibration analyses of built-up structures at high frequency ranges have received much attention because of technological advances and the pursuit of comfort.However, at high frequency ranges, the traditional noise and vibration analytical tools such as the finite element method (FEM), the boundary element method (BEM) and so on, are ineffective because of the requirement that the mesh needs to be very fine which makes model generation difficult, turn-around time longer, and computational cost relatively high.Additionally, because of increasing uncertainty of various analytic parameters in structural models at high frequencies, the confidence of analytic results obtained by the traditional tools decreases.Statistical energy analysis (SEA) is one of the methods that were developed to overcome the weaknesses of the traditional vibration tools at high frequency ranges.SEA predicts modal-and frequency-average vibrational behavior of a structure and uses energy as a primary variable [13].However, SEA is only well suited for high frequencies, where the modal density is high, and cannot predict the spatial variation of energy in a subsystem because the subsystem has a single energy value.Energy flow analysis (EFA) was developed to improve existing NVH analysis tools such as FEM, BEM and SEA.EFA models the flow of vibrational energy in a structure as heat flow and uses a heat conduction-type differential equation as governing equation.EFA can predict the spatial variation of energy and consider the local damping treatment in a structure.
EFA was first introduced by Belov et al. [15].Recently, Nefske and Sung implemented a finite element formulation of the vibrational energy equation [2].The developed energy flow finite element method (EFFEM) was used to predict the vibrational response of Euler-Bernoulli beams excited by a point harmonic force.Wohlever, in his study of vibrational energy flow in rods and Euler-Bernoulli beams, found that the vibrational far-field intensity is proportional to the gradient of the vibrational far-field energy density when all the quantities are time-and locally space-averaged for finite Euler-Bernoulli beams [3,4].Bouthier extended the work of Wohlever to energy flow in membranes, thin plates and acoustic cavity.Bouthier found that the vibration energy flow in two-dimensional structures and acoustic cavities is also governed by a second-order differential equation as in one-dimensional structures [7][8][9][10].Cho established the finite element formulations for the one-and two-dimensional energy equations by developing the joint element matrix equation for various structural joints [11,12].Park formulated the improved vibrational analysis for complex structures by developing the energy flow model of in-plane waves in finite coupled thin plates [1].
Until now, most researches on EFA have been restricted to the analysis of simple structures such as Euler-Bernoulli beams and Kirchhoff plates in structural elements.Though EFA is a more suitable method at high than low frequency range, the traditional energy flow models were not able to consider the effects of rotatory inertia and shear distortion, which are important at high frequencies.Therefore, energy flow models for superior structural elements such as Timoshenko beams and Mindlin plates were needed to improve the EFA vibrational predictions of a structure in the high frequency ranges and to incorporate the effects of rotatory inertia and shear distortion.In relation to this work, Ichchou et al. suggested a multi-propagative, simplified energy model of the Timoshenko beam [6].However, the suggested energy governing equations were not derived in detail and were not sufficiently verified in their works.
In this paper, the energy governing equations for the corrected flexural wave, accounting for the effects of rotatory inertia and shear distortion in the Timoshenko beam, are derived in detail.According to the characteristics of waves, the propagating corrected flexural waves existing in the Timoshenko beam are classified into bending dominant flexural wave (BDFW) and shear dominant flexural wave (SDFW).Unlike the traditional energy flow models, the derived energy governing equations have different forms depending on frequency.When the excitation frequency is lower than the critical frequency determined by material properties and dimensions of the beam, the energy governing equation for only BDFW like the Euler-Bernoulli beam is obtained.When the excitation frequency is higher than the critical frequency, two energy governing equations for BDFW and SDFW are obtained.Finally, to verify the accuracy and validity of the developed energy flow model, the energy flow solutions and classical solutions for the simple finite Timoshenko beams are compared for several different conditions, and the results of the Timoshenko beam model are also compared with those of the Euler-Bernoulli beam model.

Timoshenko beam theory
The Timoshenko beam theory includes the effects of rotatory inertia and shear distortion in the Euler-Bernoulli beam model.These effects become increasingly important for shorter wavelength of a vibrational wave in the medium-to-high frequency ranges.Usually, the effect of shear distortion is more significant than that of rotatory inertia on flexural vibration [5].
The equation of the corrected flexural motion in a Timoshenko beam can be obtained by using the Hamilton's variational principle.The Lagrangian is defined by the subtracting the potential energy from the kinetic energy.In a Timoshenko beam, the kinetic energy can be classified into translational kinetic energy and rotational kinetic energy; where v is the transverse displacement, α is the angle of rotation due to bending, A is the cross sectional area of the beam, I is the second moment of area of the beam's cross section, and ρ is the density of the beam.The potential energy is classified into potential energies associated with bending moment and shear force; where E is the Young's modulus, κ is the shear factor, and G is the shear modulus.
Using Eqs (1), ( 2), ( 3) and ( 4), the Lagrangian is given by Using the Hamilton's variational principle in Eq. ( 5), the equations of corrected flexural motion and the boundary conditions are obtained by, respectively,

Group velocity of corrected flexural wave in the Timoshenko beam
The general solutions for the homogeneous undamped problem of corrected flexural motion in a Timoshenko beam can be expressed by where C is the constant coefficient, d is a vector of constant numbers corresponding to an eigenvector, K is the wavenumber of a corrected flexural wave and ω is the excitation frequency.When Eq. ( 8) is substituted into Eq.( 6), the following matrix equation is obtained; In order to have a non-trivial solution, the following characteristic equation is obtained in the upper matrix equation; By Eq. ( 10), the four roots are given by Of the four roots, two are always imaginary and the other two roots, are either real or imaginary depending on the frequency ω.The frequency determining whether these roots are real or imaginary is called the critical frequency ω c and it can be expressed by For a steel beam (B × H = 0.01 m × 0.01 m), which has a rectangular cross-section, the critical frequency is about 15.8 kHz.
The corresponding eigenvectors d i can be expressed by In Eq. ( 11), the term of ρω 2 A EI, which is equivalent to the fourth power of the bending wavenumber in the Euler-Bernoulli beam can be regarded as that related to the effect of bending and the term of can be regarded as that related to the effect of shear distortion.When ω ω c , the term having the bending effect (ρω 2 A EI) becomes dominant to that having the shear effect (ρ 2 ω 4 (1/E − 1/κG) 2 4) in Eq. ( 12), as the frequency decreases.When ω > ω c , as shown in Fig. 1, the term having the shear effect (ρ becomes dominant to the term having the bending effect (ρω 2 A EI) in Eqs ( 12) and ( 13), as the frequency increases.Therefore, below the critical frequency (ω ω c ), only the wave with wavenumber k 1 propagates and its wavenumber can be approximated as where K 1 = ±k 1 j.The wavenumbers k 1 and k 2 of the propagating waves above the critical frequency (ω > ω c ) can be approximated as, respectively, where K 1 = ±k 1 j and K 2 = ±k 2 j.
In Eqs ( 16) and ( 17), the flexural wave with wavenumber k 1 is influenced by bending below the critical frequency and by shear above the critical frequency.On the other hand, the flexural wave with wavenumber k 2 propagates in only the upper frequency ranges where the effect of shear is dominant, than the critical frequency.By these characteristics of wavenumbers, the flexural waves with wavenumbers k 1 and k 2 can be called bending dominant flexural wave (BDFW) and shear dominant flexural wave (SDFW), respectively.
Therefore, the two kinds of group velocities of propagating flexural waves can be obtained.Firstly, when ω ω c , one group velocity is obtained because one kind of propagating flexural wave exists.Secondly, when ω > ω c , two group velocities are obtained because two kinds of propagating waves exist.
The group velocity of BDFW can be obtained by the derivative Similarly, the group velocity of SDFW is obtained; In case of the Euler-Bernoulli beam, the group velocity is given by where c gE is the group velocity for a flexural wave and k f is the flexural wavenumber ( ω 2 ρA EI 1/4 ) in an Euler-Bernoulli beam.
Figure 2 shows the values of group velocities in the Euler-Bernoulli beam and the Timoshenko beam with variation of frequency.

Above the critical frequency
When ω > ω c , two kinds of corrected flexural waves, that is, four propagating waves exist.These waves are commonly referred to as the far-field solution.The exponentially decaying evanescent waves, referred to as the near-field solution, can be safely neglected in the solution, though they do not exist above the critical frequency [4].Using one of the eigenvectors, which are given by Eq. ( 15), the general solutions with hysteretic damping can be expressed by where G c and E c are the complex shear and Young's modulus, which are defined as G c = G (1 + jη) and E c = E (1 + jη), respectively, and η is the structural (hysteretic) damping factor, The far-field total energy density is the sum of the potential and kinetic energy densities.The time averaged far-field energy density can be expressed by where operator indicates a time averaged quantity and the asterisk notation indicates a complex conjugate.The far-field total power in the corrected flexural motion is transmitted by the shear force and the moment.The power associated with the shear force (q s ) is where V y is the shear force (V y = κG c A (∂v/∂x − α)).
The power associated with the moment (q m ) is where M y is the bending moment (M y = E c I (∂α/∂x)).
The time averaged far-field power can be expressed by In lightly damped structures, i.e., η << 1, using Eq. ( 17), the real and imaginary components of the complex wavenumber k 1c = k 11 + jk 12 , k 11 , and k 12 , respectively, are well approximated as For the other complex wavenumber k 2c , using Eq. ( 18), the real and imaginary components of k 2c = k 21 + jk 22 , k 21 and k 22 , respectively, are well approximated as Substituting Eq. ( 22) into Eqs (24) and ( 27), the time averaged far-field energy density and power of the corrected flexural wave are obtained.Assuming that the structural damping for corrected flexural motion of the beam is sufficiently small, the terms on the order of square or more of η can be neglected.Additionally, the locally spaced average is used to eliminate the terms which are spatially harmonic in the time averaged far-field energy density (Eq.( 24)) and power (Eq.( 27)).By this procedure, the time-and locally space-averaged far-field power and energy density of the corrected flexural wave in the Timoshenko beam can be approximated as (Appendix A) where means the time-and locally space-averaged quantity.As shown in Eqs (32) and (33), above the critical frequency, the far-field power and energy density of the corrected flexural wave in the Timoshenko beam are separated into the terms of each wave component, which is BDFW or SDFW.

Below the critical frequency
When ω ω c , only two far-field solutions exist and the other solutions are near-field solutions.The far-field solutions neglecting the near-field solutions, are given by where In lightly damped structures, i.e., η 1, using Eq. ( 16), the real and imaginary components of the complex wavenumber k 1c = k 11 + jk 12 , k 11 , and k 12 , respectively, are well approximated as As above the critical frequency, by substituting Eq. (34) into Eqs ( 24) and ( 27), the time-averaged far-field energy density and power of the corrected flexural wave are obtained.Assuming that the structural damping is sufficiently small, the terms on the order of the square or more of η are neglected.The locally spaced average is also used to eliminate terms that are spatially harmonic in the time averaged far-field power and energy density.The timeand locally space-averaged far-field power and energy density in the Timoshenko beam can be approximated as (Appendix A)

Above the critical frequency
For a lightly damped rod, a one-dimensional structure, a simple relationship exists between the local values of power and the gradient of the energy density.In case of the Euler-Bernoulli beam, the relationship between the local values of far-field power and energy density of a flexural wave could be obtained by the time average and local space average [4].No such relationship between the local values of far-field power and energy densities of the corrected flexural wave in the Timoshenko beam can be easily developed for the general case, unlike the Euler-Bernoulli beam.
When ω > ω c , the time-and locally space-averaged far-field power and energy density, which are expressed by Eqs (32) and (33), respectively, are composed of two kinds of wave components, the BDFW and SDFW.Both far-field power and energy density can be separated into the terms of each wave component.The time-and locally space-averaged far-field energy density can be expressed by ē = ē1 + ē2 and (40) The time-and locally space-averaged far-field power can be expressed by In Eqs (40) and (41), "1" and "2" in the subscripts denote BDFW and SDFW respectively.The time-and locally space-averaged far-field energy density and power of BDFW with wavenumber k 1c are given respectively by Above the critical frequency, because the effect of shear is dominant to that of bending, the square of k 11 , which is the real part of the wavenumber k 1c , can be approximated as Using Eq. ( 46), the group velocity (Eq.( 19)) of BDFW having the wavenumber k 1c can be approximated by Substituting Eq. (46) into Eqs (42) and (43), the relationship between the time-and locally space-averaged far-field power q1 and energy density ē1 , is derived by In case of SDFW having the wavenumber k 2c , by the procedure used in the previous case, the relationship between the time-and locally space-averaged far-field power q 2 and energy density ē2 can be obtained.Above the critical frequency, because the effect of shear are dominant to that of bending, the square of k 21 can be approximated as Using Eq. ( 49), the group velocity (Eq.( 20)) of SDFW having the wavenumber k 2c can be expressed by Substituting Eq. (49) into Eqs (44) and ( 45), the relationship between the time-and locally space-averaged far-field power q2 and energy density ē2 , is derived by The power injected by the external loads into an elastic medium is dissipated by damping and is transmitted to the neighboring media.For a steady state elastic system, the power balance equation can be written as where π diss and π in are the dissipated power by the damping of the system and the input power, respectively.From the works of Cremer and Heckl [5], the time averaged dissipated power in an elastic medium with small structural damping is proportional to the time averaged total energy density in the form By using the power balance equation (Eq.( 52)), the energy transmission relations (Eqs (48) and ( 51)) and the power dissipation relation (Eq.( 53)), the second order partial differential equations which take the far-field total energy density related to each kind of wave as a primary variable, can be derived respectively, The homogeneous solutions of the two energy governing equations are obtained respectively as follows.
ē1 = Ae −φ1x + Be φ1x and (56) where A, B, C and D are constants, which are determined by the boundary conditions, and φ 1 and φ 2 are defined as ηω/c g1 and ηω/c g2 , respectively.Considering the global energy density (Eq.( 40)) of the sum of ē 1 and ē2 , the global energy density can be expressed by To express the solution as Eq. ( 58), a fourth order differential equation that takes the global energy density as a primary variable is required.This homogeneous differential equation can be expressed as Multiplying the left side of Eq. ( 59) by (c g1 c g2 ) 2 (ηω) 3 to express the quantity of each term in Eq. ( 59) as that of a power, the following new-typed energy governing equation can be obtained as Using Eqs (52) and (60), the relationship between the time-and locally space-averaged global energy density and global power can be expressed as Naturally, this energy transmission relation (Eq.( 61)) can be verified by substituting Eqs (32), ( 33), ( 46), ( 47), ( 49) and (50) into Eq.( 61).The energy governing equation which is analogous to Eq. ( 60) was simply suggested by Ichchou et al. [6].However, they did not directly derive the energy transmission relation (Eq.( 61)), which is most important in the derivation of the energy governing equation, for the corrected flexural wave in the Timoshenko beam by using the classical solutions based on the displacement solutions for the Timoshenko beam and did not sufficiently verify suggested energy governing equations.In their works, assuming that the energy transmission relations (Eqs (48) and ( 51)) of the two kinds of flexural waves in the Timoshenko beam are established, the relationship between global energy density and global power (Eq.( 61)) was obtained by simple combination of the two assumed energy transmission relations.However, though the energy transmission relations in the rod and Euler-Bernoulli beam are found [4], the establishment of the energy transmission relations (Eqs ( 48) and ( 51)) in the Timoshenko beam can not be guaranteed without direct derivation in Section 5. Additionally, the fourth order differential energy governing equation (Eq.( 60)), which needs four boundary conditions in order to obtain its solution, is not an appropriate form for describing the energetics of the corrected flexural wave in the Timoshenko beam.In one-dimensional energetic models such as those of the Timoshenko beam, the model having a homogeneous domain has only two boundary conditions at both ends.Therefore, for the corrected flexural wave of the Timoshenko beam, the second order differential energy governing equations (Eqs (54) and ( 55)) representing the energetics of each kind of propagating wave are more reasonable than the fourth order differential energy governing equation (Eq.( 60)) representing the global energy of two kinds of propagating waves.Above all, to analyze the coupled Timoshenko beam structure as will be presented later in the other companion paper, the energetics of one kind of propagating wave like Eq. ( 54) or Eq. ( 55) has to be considered.

Below the critical frequency
When ω ω c , one kind of BDFW having the wavenumber k 1c exists.The time-and locally space-averaged far-field energy density and power are given by Eqs (39) and ( 38) respectively.Below the critical frequency, because the bending effect is dominant than shear effect, the square of k 11 , which is the real part of k 1c , can be approximated as Using Eq. ( 62), the group velocity (Eq.( 19)) of BDFW having the wavenumber k 1c can be expressed by This approximated group velocity is identical to that of the flexural wave in the Euler-Bernoulli beam.Substituting Eq. (62) into Eqs (38) and (39), the relationship between the time-and locally space-averaged far-field power q and energy density ē , is derived by By using the power balance equation (Eq.( 52)), the energy transmission relations (Eq.( 64)) and the power dissipation relation (Eq.( 53)), the second order partial differential equation which takes the total energy density as a primary variable can be derived as When ω ω c , the energy governing equation for a BDFW in the Timoshenko beam is almost equivalent to that for a flexural wave in the Euler-Bernoulli beam.Therefore, the derived energy governing equations are represented in the general form including not only the Euler-Bernoulli beam theory used for the conventional energy flow model [4] but also the Rayleigh, shear, and Timoshenko beam theories.

Numerical examples
To verify the derived energy governing equation for the corrected flexural wave of the Timoshenko beam, numerical analyses are performed for the finite Timoshenko beam simply supported at both ends and excited by a transverse harmonic point force, as shown Fig. 3.The time-and locally space-averaged far-field energy density of corrected flexural wave is calculated by solution (Eq. ( 66)) of the energy governing equation developed for the vibrating Timoshenko beam, and the power distribution in the beam is obtained by Eq. ( 67).

ē = (Ae
x + Be ηω c g1 x ) + (Ce For some examples, the energy density and power distributions obtained by the developed energy flow model for the Timoshenko beam are compared with the results obtained by the classical solution for the Timoshenko beam, the classical solution and the EFA solution for the Euler-Bernoulli beam [4].In the first example, the dimensions of the beam, shown in Fig. 3, are L × B × H = 1 m × 0.01 m × 0.01 m and the material properties of the beam are assumed to be the same as those of steel.The force is located at x 0 = L/2 in the beam and its amplitude is F = 1 N.The time averaged input power for the EFA solutions can be calculated as follows: where ∂w/∂t can be obtained by the classical solution expressed as Eq.(B4) or (B5).The slenderness ratio (L/H) of this model is 100, and the critical frequency is about 159 kHz.When the excitation frequency is f = 1 kHz, the spatial distributions of the energy density and power obtained by each solution for various values of structural damping (η = 0.01, η = 0.05, η = 0.1 and η = 0.2) are shown in Figs 4 and 5, respectively.Figures 4 and 5 show that the developed EFA solutions represent well the global variation of the classical solutions for the Timoshenko beam, regardless of structural damping.The beam model in this example can be regarded as a thin beam, because of the large value for the slenderness ratio (L/H = 100) and the high critical frequency (f c = 159 kHz).When f = 1 kHz, because the wavelength of the BDFW (λ ≈ 0.3) is longer than the maximum length of the cross-section of the beam (H = B = 0.01), the effects of rotatory inertia and shear distortion are not dominant to the effect of bending.Therefore, as shown in Figs 4 and 5, the results obtained by the Timoshenko beam theory are almost equivalent to those by the Euler-Bernoulli beam theory.Figures 6 and 7 show the distributions of the energy density and power obtained by each solution for various excitation frequencies when the structural damping is η = 0.01.In Figs 6 and 7, the developed EFA solutions of this case, as in the previous case, represent well the global variation of the classical solutions for the Timoshenko beam, regardless of the excitation frequency.On the other hand, as the excitation frequency increases, the results from the Timoshenko beam model and the results from the Euler-Bernoulli beam model become widely different because the effects of rotatory inertia and shear distortion, which are

Conclusion
For the improved vibrational analysis of beam structures in the medium-to-high frequency ranges, the energy flow model for the finite coupled Timoshenko beam was newly developed.First, the energy governing equations  Rayleigh, shear, and Timoshenko beam theories.When the excitation frequency is lower than the critical frequency, the energy governing equation for one propagating flexural wave, BDFW, like the energy governing equation of the Euler-Bernoulli beam is obtained.When the excitation frequency is higher than the critical frequency, the energy the excitation frequency and damping loss factor increase, the results from the energy flow models between the Timoshenko beam and Euler-Bernoulli beam were remarkably different.
The derived energy flow model for the finite Timoshenko beam can be effectively used for the advanced prediction of the time-and locally space-averaged far-field energy density and power distributions of built-up structures composed of beams in the high frequency ranges if the extended researches on the coupling relationships of the Timoshenko beam are performed, which are described in the other companion paper.
Using the classical solutions solved by these boundary conditions, the time averaged total energy density and power are obtained by Eqs (24) and (27).

Fig. 9 ,
Fig.9, the effects of rotatory inertia and shear distortion are dominant in almost all frequency ranges.Therefore, as the excitation frequency increases and the beam grows thicker, the effects of the rotatory inertia and shear distortion in the Timoshenko beam become very important.