High-Frequency Dynamic Analysis of Plates in Thermal Environments Based on Energy Finite Element Method

The energy density governing equation to analyze the high-frequency dynamic behavior of plates in thermal environments is derived in this paper, in which the thermal effects are considered to change the membrane stress state and temperature dependent material properties of plates. Then the thermal effects on the energy reflection and transmission coefficients are dealt with hereof. Based on the above, an EFEM (energy finite element method) based approximate approach for the energy analysis of coupled plates under nonuniform thermal environments is proposed. The approach could be conducted by three steps: (1) thermal analysis, (2) thermal stress analysis, and (3) forming element matrixes, joint matrixes, and the whole EFEM formulation for the energy analysis. The same mesh model is used for all the three steps. The comparison between EFEM results and classical modal superposition method results of simply supported plates in various uniform thermal environments and coupled plates in nonuniform thermal environments demonstrated that the derived energy governing equation and the proposed approach described well the smooth time- and locally space-averaged energy density. It is found that the distributions and levels of energy density are affected by thermal effects, and the variation trends are related to exciting frequency.


Introduction
With the great development of hypersonic crafts which are usually subjected to extremely aerodynamic heating and high-frequency exciting during working, there is a great need for the high-frequency dynamic analysis of structures in thermal environments.
Thermal environments have a series of effects on material properties, geometry shapes, stress state, and so on.Great efforts have been made to analyze the dynamic characters of the structure in thermal environments at low frequencies.Ganesan and Dhotarad [1] developed a numerical method for the vibration analysis of thermally stressed plates in which the thermal stresses were evaluated by the finite element method and these stress values were then used in a dynamic analysis of the plate performed by either the finite difference method or variation methods.Jeyaraj et al. [2,3] used combined FEM/BEM to analyze the vibration and acoustic radiation characters of an isotropic plate and a composite plate with inherent material damping in a thermal environment, finding out that the natural frequencies decrease when the temperature increases and the overall sound radiation of the plate reduces marginally only due to the interaction between reduced stiffness and enhanced damping.Geng and Li [4] analyzed the acoustic radiation and vibration of a flat plate in thermal environments through both the theory method and combined FEM/BEM, and the two results correspond to each other.Liu and Li [5] investigated the vibration and acoustic response of rectangular sandwich plate in thermal environments.
With the frequency increase, more and more elements are needed to describe the vibration of structures, and also small uncertainties will have more and more effects on the calculating results.Thus the FEM and BEM calculation will not only cost much time, but the error is much bigger.Though SEA (statistical energy analysis) is widely used to predict the space-and frequency-averaged behavior of builtup structures at high frequencies where the modal density of 2 Shock and Vibration structures is high [6], it could only approximate a single acoustic or vibrational energy value for each subsystem of structures.
Due to the disadvantages of the FEM, BEM, and SEA at high frequency, EFEM (energy finite element method) which could predict the time-and space-averaged far field vibrational energy of structures appeared.Belov et al. [7] first proposed the power flow model.Nefske and Sung [8] developed the finite element formulation of the power flow model and applied it to beams.Wohlever and Bernhard [9] made a further research on the vibration energy response of rods and beams.Bouthier and Bernhard extended the work to two-dimensional structures and analyzed the energy density distribution of a membrane [10] and a Kirchhoff plate [11].Wang et al. [12] presented the high-frequency energy boundary element method combined with energy finite element to describe the energy density distribution of structure-acoustic coupled system.Park et al. [13] developed the power flow model of in-plane waves in thin plates and flexural waves in finite orthotropic plates.Then they derived the energy governing equation of flexural waves in Timoshenko beam [14,15], Mindlin plate [16], and Rayleigh-Bishop rod [17] which take shear distortion and rotatory inertia acting an important role in high-frequency range into consideration.Zhang et al. developed an alternative energy finite element formulation for interior acoustic spaces and thin plates considering the wave response as a summation of incoherent orthogonal waves [18].They also expanded the EFEM to analyze the energy distribution of stiffened plates under heavy fluid loading [19].Xie et al. applied EFEM to high-frequency structural-acoustic coupling of an aircraft cabin with truncated conical shape [20] and proposed the transient vibrational energy response analysis of a rod under high-frequency excitation [21].Park made an introduction of the developed EFEM based software-EFADS, and the results of one part of a real automobile derived from EFADS matched well with the experimental results [22].
Though EFEM has been developed well since 1980s, only Zhang et al. [23] researched on the thermal effects on the high-frequency vibration.They derived the energy density governing equation of beams and verified the accuracy.In this paper, we derived the energy governing equation of plates in thermal environments.The thermal environments are considered to affect the material property and thermal stress condition.We then study the thermal effects on power transmission and reflection coefficient and develop an approximate approach to analyze the energy density distribution of coupled plates in nonuniform thermal environments: first, the thermal analysis and thermal stress analysis are conducted; then based on the temperature and thermal stresses derived above, the energy distribution of the structure could be obtained therefore.The accuracy of the derived equation is verified by comparing the developed EFEM results with classical modal superposition results of a simply supported plate in uniform thermal environments.The thermal effects on the energy distribution and level are also analyzed.Finally, the numerical example of the coupled plate in nonuniform thermal environments shows that the approach could describe well the thermal effects on the level and distribution of energy density.

Energy Governing Equation for a Plate in Thermal Environments
2.1.Derivation of Wavenumber and Group Velocity with Thermal Effects.For a plate with the temperature uniform throughout the thickness, assume that the plate is stressfree at the reference temperature  0 .With the temperature increasing, the changes of thermal stress and material properties should be taken into consideration.So when the temperature changes to , the transversely vibrating governing equation will be [4] where  = ℎ 3 (1+)/12(1− 2 ) =  0 (1+) is the complex bending stiffness of the plate,  is the structural damping loss factor,  is the density, ℎ is the thickness,  is the transverse mechanical load applied on plate, and  is the imaginary unit.  ,   , and   are the thermally induced membrane forces.
Use the separation of variables, the general solution can be expressed as where   and   are the wavenumbers in  and  directions and  is the circular frequency.At small damping, the complex wavenumbers in  and  directions above could be approximate as By substituting (2) into (1), the dispersion relation can be derived as Assume that  is the total wave numbers in the wave propagation angle , so  1 and  1 could be expressed as [25]  1 =  cos   1 =  sin .
Substituting ( 5) into (4) Assume that Thus we could derive the expression of  in terms of : From ( 8) and (7) we could find that the total wavenumber  varies with the wave propagation angle .In EFEA method, the plate is diffuse wave fields, so averaging  1 and  1 over , we could obtain the averaged diffuse wavenumber in  and  directions: The group velocity in terms of  could be expressed as Equation (10) shows that the group velocity also depends on the direction of the wave propagation.So the same as ( 9), the group velocity should also be averaged over  [26]: From ( 9) and (11), we can see that the wavenumber and group velocity will be changed in thermal environments.Because, in high-frequency range, the wavelength is quite small and the near field waves usually decay in half wavelength, we only utilize the far field solution of (1) to simplify the energy analysis [11].So the general form of the far field travelling plane wave solution can be expressed as where the unknown constants   ,   ,   , and   are the amplitudes of the waves in  and  directions.

Derivation of the Energy Governing Equation with Thermal Effects.
As discussed by Bouthier and Bernhard [11], the time-averaged energy density which includes both the kinetic and potential energy densities of the plate across the thickness can be expressed in terms of transverse displacement as As above, the time-averaged energy intensity across the thickness can be expressed in terms of transverse displacement as [11] Substituting ( 12) into ( 13)-( 14), we can obtain the complete time-averaged far field energy density and energy intensity expressions in terms of transverse displacement.For we could not find obvious relations between them, they are averaged spatially over a half wavelength for small damping in the following manner [11]: Thus the time-and space-averaged density ⟨⟩ can be written as The time-and space-averaged intensities ⟨  ⟩ and ⟨  ⟩ can be written as Obviously the relationship between the time-and spaceaveraged energy density and intensity is For an elastic medium, the differential steady energy balance equation in terms of time-and space-averaged energy density can be written as [8] Combining ( 18) and ( 19), we can obtain the steady far field energy governing equation of plates in uniform thermal environments: It is obvious that thermal environments affect the first coefficient of (20) by changing wavenumbers of the plate; thus the energy distribution state of the plate is changed thereof.Compared to the energy governing equation derived by Bouthier and Bernhard [11], (20) becomes the same when the temperature becomes reference temperature.
The finite element formulation is introduced here to solve (20) [19]: where   is the system matrix for element including thermal effects,   is the vector of nodal energy density for the corresponding element,   is the vector of input power at the nodal locations, and   is the vector of power flow across element boundaries including thermal effects, which will be introduced in the next sections.
In high-frequency range, for the wavelength is small, the input power could be expressed as [27] or where  is the amplitude of the external point force,  * is the conjugate of the velocity of the exciting point and  is the impedance of the corresponding infinite plates with thermal effects at the driving point expressed as follows: The wavenumber will be anisotropic when the thermal stresses are different in different directions, and the structure is diffused field, so the wavenumber is averaged over .

An Approximate Method for the Energy Density Analysis of Coupled Plates in Arbitrary Thermal Environments
Usually, the thermal environments are nonuniform, and the structures are made up of coupled plates, which means the thermal stresses and temperature dependent material properties are not the same everywhere.Thus the group velocities will be different everywhere.And also the power transmission and reflection will happen everywhere.Therefore the derived energy governing equation ( 20) cannot be used directly for structures as traditional method.So it is of great importance to develop a method to analyze the high-frequency vibration energy density of coupled plates in arbitrary thermal environment.

Derivation of Power Transmission and Reflection Coefficients in Thermal
Environments.As we know, it is not easy and takes a lot of time to consider all the energy transmission and reflection.So we assume that the temperature and thermal stresses are uniform in every element and just calculate the energy transmission and reflection coefficients between every two adjacent elements.The method introduced by Langley and Heron [28] is expended here to take the thermal effects into consideration for the plates.As shown in Figure 1(a), the model of plate junction in global coordinate is introduced here.The tractions (  ,   ,   ), displacements (  , V  ,   ), thermal stresses (  ,   ,   ), and temperature   of either single plate are defined in local coordinate as shown in Figure 1(b).For the membrane stress will not affect the in-plane vibration characteristics, the thermal effects will change the in-plane energy density distribution and level only by the temperature dependent material properties.So the governing equations of out-of-plane and in-plane deformation of the th plate defined in local coordinate are z j (S j w j ) x j (T j u j ) The relationship between the forces and moments per unit length which are applied to the junction by plates and the tractions which act at the edge of plates could be expressed as where the forces and moments vector Q are defined in global coordinate and the tractions vector The compatibility conditions between the junction displacements a = ( V  ) and the edge displacements of plate  b  = (  V      )  could be expressed as According to (25), the wavenumbers corresponding to bending, longitudinal, and shear waves are Assume that the incident wave travels through one of the plates towards the junction with an incident angle of   .The compatibility at the junction means the incident waves along with the reflected waves and transmitted waves share the same wave motion  − 1 sin    in  direction.So the wave motion in  direction could be obtained by where the subscript  means the number of an arbitrary plate.So the bending and in-plane response could be expressed as where   ,   , and   are the complex amplitudes of each wave.
According to [28], the equation determining the junction displacements is where K  is the dynamic stiffness matrix of th plate.The th plate carries the incident wave.By solving (32), the displacement of junction in general coordinate is derived.And the amplitudes of each wave could be obtained through (28) and (31).And only propagating waves carry power away from the junction.The power could generally be written as 2  sin , where  is the average kinetic energy and   is the group velocity.So the power transmitted by each type of wave could be expressed as where  is the propagating angle of corresponding wave.The expression of   is different from that in [28] for the group velocity has been changed in thermal environments.
The power transmission coefficients    (, ) associated with each of the generated waves are defined as the ratio of power of generated waves over power of incident wave, where , , , and  represent the carrier plate, wave type, frequency, and heading of the incident wave and  and  represent the carrier plate and wave type of a generated wave.The diffuse field transmission coefficients are calculated by Therefore, the joint matrix in EFEA formulation is [29] where   and   are Lagrangian basis functions and  is the boundary between element  and element .Substituting (35) into ( 21), the whole EFEM formulation is derived and the thermal effects are taken into consideration completely.

An Approximate Analysis Approach to Analyse the Energy Distribution of Plates in Nonuniform Thermal Environments.
Based on the energy governing equation (20) and the reflection, transmission coefficients in thermal environments derived above, we develop the numerical approach to analyze the energy distribution of coupled plates in arbitrary thermal environments.
Thermal analysis and thermal stress analysis only need a few elements and the same mesh model could also be used for energy density distribution analysis.So as shown in Figure 2, in the complete analysis process, thermal analysis and thermal stress analysis should be done first, and the four nodal temperatures and thermal stresses of every elements are averaged to calculate the reflection and transmission coefficients, element matrixes, and joint matrixes.Then the whole EFEM matrix is formed to calculate energy density distribution of the whole structure.

Validation of the Derived Energy Governing Equation for Simply Supported Plates in Uniform Thermal Environments.
In this section, we will verify the accuracy of the new derived energy governing equation of simply supported plates in uniform thermal environments and investigate how the thermal environments affect the energy distribution of the plate.
As shown in Figure 3, a simply supported plate is studied here.The plate is excited by a harmonic point force at the center.The geometric parameters and temperature independent material properties are shown in Table 1.The temperature dependent elastic modulus and thermal expansion coefficient especially are presented in Table 2 [24].
According to the boundary condition and material properties above, the thermal stress could be calculated as That means the wavenumbers and the group velocities are all the same in different directions.Assuming that   =   = , the wavenumber and group velocity could be simplified as Thus the energy governing equation of the simply supported plate in uniform thermal environment could be simplified as    20) is developed here to calculate the energy density distribution of the example.The results are compared to the classical modal superposition results.As discussed by Geng and Li [4], the displacement modal superposition expression of (1) with four boundaries simply supported and a harmonic point force excited in the center is where Thus the exact energy density could be obtained by substituting (39) into (13).
The input power is computed from (22), and the velocity could be obtained from the modal superposition method.
First we choose 4 different cases here to verify the accuracy of the derived energy governing equation.The temperature, frequency, and damping factor of each case are shown in Table 3.The reference temperature in this work is 20 ∘ C.And the reference energy density is 1 × 10 −12 J/m 2 .Figures 4, 5, 6, and 7 give the calculated energy density distributions of the four cases obtained by both the classical modal superposition method and EFEM, respectively.It is clear that the classical modal solutions are fluctuating surfaces, while the EFEM solutions are smooth surfaces which represent the averaged energy density distribution trends.And from the observation of the four figures, we can find that the global variation of energy density increases with the damping factor increasing, and the energy density level decreases with the frequency increasing.
The comparisons of the energy density along the line  =  of the four cases are illustrated in Figures 8 and 9.It is obvious that in all the four cases, EFEM solutions approximate well the time-and space-averaged far field energy density compared to the classical modal solutions.
From the above cases, the good agreement between EFEM solutions and modal superposition solutions demonstrates that the energy governing equation derived in this paper approximates well the smooth time-and space-averaged energy density distribution of simply supported plate in thermal environments.

Investigation of Thermal Effect on Energy Density Level.
In this section, the effects of thermal environments on the averaged energy density level of the plate are investigated.The energy density levels obtained by EFEM in 961 points uniformly distributed in the plate are averaged as the averaged energy density level of the plate.
First, the averaged energy level variation by exciting frequency from 5000 Hz to 50000 Hz with an increment of 100 Hz when temperatures are reference temperature 20 ∘ C and 65 ∘ C is shown in Figure 10.As the figure shows, both in 20 ∘ C and 65 ∘ C, the energy density level curves decline with the frequency increasing on the whole, while fluctuating in every local region, and at the same time, the curve in 65 ∘ C wholly moves to the left compared to the curve in 20 ∘ C.   The reason why both curves decline with fluctuating is that the vibration energy of the structure is mainly distributed at lower frequency range and will rise suddenly at natural frequencies for resonance while the density of natural frequencies increases with the frequency increase.As for the movement of curve in 65 ∘ C, it could be explained as the thermal effects: the thermal stress and the decrease of elastic modulus both result in the reduction of bending stiffness which make the natural frequencies wholly reduce at the same time.Then the local energy density level peaks will wholly float to lower frequencies.
Above all, the variation trend of energy density according to temperature depends on the exciting frequency and the temperature increasing range.For the case of a small increasing range, the frequency response curve only moves a little that no peaks or valleys pass the exciting frequency; then the energy density will increase or decrease monotonously according to the exciting frequency.For the case of a big increasing range, which means that the curve moves a lot that any peaks or valleys pass the exciting frequency, then the variation trend of the energy density will be nonmonotonous.

Numerical Examples of Coupled Plates in Nonuniform
Thermal Environments.The model here is the simply supported coupled plates as shown in Figure 11, where the joint angle is right angle and the dimension and material properties   of the two plates are the same as that in Tables 1 and 2. Table 4 gives the detail parameters of the two cases.The temperature boundary condition of case 1 especially is that the temperatures of the centers of the two plates are set to 80 ∘ , while the temperatures of the boundaries of the coupled plates are the reference 20 ∘ .Case 2 is set to the reference uniform thermal environment as a comparison with case 1.
The transverse harmonic point force is located at the center of plate 1 with an amplitude of 10 N and an exciting frequency of 15000 Hz.
For case 1, the plate is divided into 400 uniform elements.Thermal analysis and thermal stress analysis are carried out first to obtain the temperature field and thermal stress field of the plate.The input power is computed from (22), and the velocity could be obtained from FEA dynamic analysis.On the base of the temperature field and thermal stress field, the   FEM solution is illustrated in Figure 12.Partially enlarging Figure 12, energy jumps caused by thermal effects will be observed at the joints of every two adjacent elements as shown in Figure 13.So above all, it is obviously observed that the method proposed in this paper could not only and membrane stresses of the plate.Thus the thermal effects are included in the derivation of energy governing equation through wavenumber and group velocity which could be expressed as functions in terms of material properties and membrane stress.Then the thermal effects on the power transmission and reflection coefficients are also investigated.
Based on the governing equation and power transmission coefficients derived above, an approximate approach to analyze the energy density distribution of coupled plates in nonuniform thermal environments is developed, in which, thermal analysis and thermal stress analysis are completed firstly; then the obtained temperature and thermal stress are averaged at element level to form the element matrixes, joint matrixes, and the whole EFEM formulation, and the energy distribution could be calculated therefore.
To verify the accuracy of the derived governing equation, a series of numerical examples of simply supported plates in uniform environments are performed with both EFEM and modal superposition method.The results demonstrate that the derived energy governing equation approximates well the time-and space-averaged energy density.The results also show that the averaged energy density of the plate will vary as the temperature changes, and the variation trends depend on the exciting frequency for sudden rises will appear at natural frequencies.For the examples of coupled plates in nonuniform thermal environments, energy jumps will occur at the joints of every two adjacent elements, and thus the variation of energy density distribution caused by thermal effects could be reflected.
Hypersonic aircrafts which are made up of a lot of coupled plates usually work in nonuniform thermal environments, subjected to high-frequency exciting.The present study could be a foundation for researchers to do further research on the high-frequency vibration energy response of complex structures in nonuniform thermal environments.

Figure 1 :
Figure 1: (a) Schematic of plate junction in thermal environments.(b) Coordinate system, displacements, tractions, temperature, and thermal stresses of plate .

Figure 4 :
Figure 4: The energy density distributions of the plate in case 1 calculated by both modal superposition method and EFEM.

Figure 5 :
Figure 5: The energy density distributions of the plate in case 2 calculated by both modal superposition method and EFEM.

Figure 6 :
Figure 6: The energy density distributions of the plate in case 3 calculated by both modal superposition method and EFEM.

Figure 7 :
Figure 7: The energy density distributions of the plate in case 4 calculated by both modal superposition method and EFEM.

Figure 8 :
Figure 8: The energy density distribution comparisons between the two methods along the line  =  in case 1 and case 2.

Figure 10 : 1 Plate 2 F𝛿
Figure 10: Averaged energy density levels by frequency with an increment of 100 Hz when  = 20 ∘ C and  = 65 ∘ C.
Energy distribution analysis process of coupled plates in nonuniform thermal environments.

Table 1 :
Geometric parameters and material properties of the plate.

Table 3 :
Parameters for each case.

Table 4 :
Parameters for each case.