Introduction

Variable thickness plates have been widely used in aerospace, ship structures, and building fields, due to their excellent mechanical performance and higher structural stability than uniform-thickness plates1. When applied to transport equipment, due to the increase of running speed, variable thickness plates are highly susceptible to be induced high-frequency vibration. Severe high-frequency vibration can easily cause structural damage. So it is necessary to study vibration characteristics of variable thickness plates.

The dynamic characteristics of variable thickness plate have been studied by many scholars. Ohga2 utilized a combination of boundary element and transfer matrix methods to solve the dynamic response of rectangular variable thickness plates. In this case, the variable thickness plate was constructed by splicing three plates with sequentially decreasing thicknesses. Some researchers3,4,5 invented a concise method for solving variable thickness plates with a variety of boundary conditions and loads through the use of an equivalent flat plate. Alam et al.6 proposed a new model for solving arched bridge-like variable thickness plates with quadratic thickness variations, which utilizes the small parameter method and the effective modulus method. Li et al.7 used the transfer matrix method to establish and solve the dynamic model of variable thickness plates, and investigated the aggregation effects of strain and kinetic energies in the resonance state of the variable thickness plate. Fu et al.8 analyzed the vibro-acoustic response of the structure by using FEM for the elastic plate (including the variable thickness portion)-cavity coupling structure, and carried out the experimental validation. Wu et al.9 analyzed the acoustic absorption characteristics of a flexible variable thickness material using experimental and FEM.

All the above studies are modal-based method to analyze the dynamic response of the variable thickness structure. However, when the structure is excited by high frequency load, the vibration modes are dense and prone to localized dynamic response. At this point, in the high-frequency range, the modal-based method faces some challenges. Take the most representative modal-based method, FEM, as an example. In FEM, the variable of governing equation is displacement, and the dynamic response is represented by modal superposition displacement. On the one side, due to the dense vibration mode, in order to obtain the structural response more accurately, smaller mesh sizes is required, which results in the more complicate computational process and the more computational time10. On the other side, due to the tendency to localized dynamic response, the structural response is extremely sensitive to design changes, which makes that FEM is not always efficient and accurate. Therefore the modal-based dynamic response analysis method is time-consuming and unable to reflect the local response information, which is not applicable to high frequencies.

The method based on energy is an effective approach for analyzing high-frequency dynamic response, typically, such as statistical energy analysis (SEA) and EFEM. In SEA, the structure is divided as several substructures. The internal response results of every substructure cannot be obtained. EFEM is a promising method to predict high-frequency dynamic responses. In EFEM, the time-and-space-averaged energy density is taken as a variable in the differential governing equation, and then the governing equations is converted into algebraic systems for solving by the Galerkin method and structure discretization. Thus energy density at all nodes in the structures can be obtained. In EFEM, the variable of governing equation is time-and-space-averaged energy density, and the dynamic response is represented by energy propagation. This is similar to the analysis of heat conduction. Therefore, when discretizing the structure, accurate dynamic response can be obtained with a smaller mesh numbers. Overall, EFEM is more suitable for high-frequency dynamic response analysis.

Nefske and Sung11 derived the EFEM governing constitutive equation for beam according to the relationship between energy flow and energy density, function expression of loss power and energy density, and relationship between power input and output in elements. By using a method similar to that of Sung, Jezequel et al.12 built the EFEM constitutive equation for longitudinal vibrating rod. Bouthier et al.13deduced the EFEM constitutive equation for thin plates and solved the equation through finite element analysis. Xie et al.14predicted the energy density distribution on cylindrical shells under axially symmetric loads based on the thin shell theory. Liu et al.15analyzed the energy dissipation relationships of beams, rods, and plates via EFEM, and deeply studied the EFEM of coupled structures.

The aforementioned researches about EFEM focus on uniform-thickness structures. In EFEM for variable thickness structures, a current approximate way (AEFEM) is to use a large number of constant thickness elements to simulate the variable thickness structures. In AEFEM, the governing constitutive equation for constant thickness is applied, which results in the inaccuracy in every element. Additionally, due to the different thickness between adjacent elements, the energy density is not continuity. So the energy transfer relationship is needed. But it will be very complicated that consider the energy transfer relationship. A compromise is that the energy transfer relationship is not considered, which will also reduces the accuracy of calculation. So AEFEM cannot accurately and completely demonstrate the characteristics of variable thickness structures and it is time-consuming16. Considering this, it is necessary to construct constitutive model for variable thickness structures and develop the NEFEM.

Xie et al.derived the NEFEM constitutive equation for variable thickness beams17 and taper beam16. It is also necessary to explore the NEFEM constitutive equation and solving method for variable thickness plates.

An NEFEM for high-frequency dynamic responses of the power-law variable thickness plate was established. Starting from the motion equation, the NEFEM constitutive equation for the variable thickness plate was constructed. In addition, the energy-density matrix equation was built using the Galerkin method and structure discretization. Through example analysis and solution, energy-density distribution of the variable thickness plate was acquired. The results was compared with the results from FEM, the NEFEM for variable thickness plate are validated. Then the advantage of the NEFEM presented in this paper is study by two comparison examples (Table 1).

The contribution to predicting of high frequency dynamic response for variable thickness plates is shown as Table 1.

Table 1 Comparison of NEFEM and current methods.
Full size table

NEFEM for power-law variable thickness plates

To study the energy finite elements method, two steps are needed. One step is deriving the constitutive equation based on energy density, and the second step is solving the constitutive equation derived.

The constitutive model for power-law variable thickness plates

The studied plate of power-law variable thickness is displayed in Fig. 1. The thickness along x direction changes according to the power law, while that along y direction is constant. The variation of thickness with x meets\(h\left( x \right)={h_0}*{\left( {1 - {x \mathord{\left/ {\vphantom {x {{x_0}}}} \right. \kern-0pt} {{x_0}}}} \right)^2}\), in which \(\:{\text{h}}_{0}\) is the thickness at the thick side of the variable thickness plate and \(\:{\text{x}}_{0}\) is the largest length of the plate along x direction. While, the length of variable thickness plates cannot reach \(\:{\text{x}}_{0}\) in practical applications, due to the machining limitation. It is assumed that the length of the variable thickness plate is \(\:{\text{x}}_{\text{a}}\).

Fig. 1
figure 1

Schematic diagram for a variable thickness plate.

Full size image

If the variable thickness plate is excited by vertically incident bending waves, the out-of-plane displacement \(w\left( {x,y} \right)\)generated due to the bending waves meets the following equation18:

$$\begin{gathered} \frac{{{\partial ^2}}}{{\partial {x^2}}}\left[ {D\left( {\frac{{{\partial ^2}w\left( {x,y,t} \right)}}{{\partial {x^2}}}+\nu \frac{{{\partial ^2}w\left( {x,y,t} \right)}}{{\partial {y^2}}}} \right)} \right]+2\left( {1 - \nu } \right)\frac{{{\partial ^2}}}{{\partial x\partial y}}\left[ {D\frac{{{\partial ^2}w\left( {x,y,t} \right)}}{{\partial x\partial y}}} \right] \hfill \\ +\frac{{{\partial ^2}}}{{\partial {y^2}}}\left[ {D\left( {\frac{{{\partial ^2}w\left( {x,y,t} \right)}}{{\partial {x^2}}}+\nu \frac{{{\partial ^2}w\left( {x,y,t} \right)}}{{\partial {y^2}}}} \right)} \right]+\rho h\left( {x,y} \right)\frac{{{\partial ^2}w\left( {x,y,t} \right)}}{{\partial {t^2}}}=0 \hfill \\ \end{gathered}$$
(1)

where \(D=E{h^3}\left( {x,y} \right)/12\left( {1 - {\nu ^2}} \right)\), \(\rho\), and \(h\left( {x,y} \right)\)separately represent the bending stiffness, density, and thickness of the variable thickness plate; E denotes the Young’s modulus of the material, v denotes the Poisson’s ratio of the material. By using the geometric acoustic approximation for one-dimensional variable thickness plates, the approximate displacement solution is expressed by

$$w\left( {x,y,t} \right)=A\left( {x,y} \right){e^{i{k_p}{S_x}\left( x \right)}}{e^{i{k_y}y}}{e^{i\omega t}}$$
(2)

where \(A\left( {x,y} \right)\)denotes the slowly changing amplitude; \({S_x}\left( x \right)\)is the part related to x in the phase function(the phase function is \(S\left( {x,y} \right)={S_x}\left( x \right)+{k_y}y\)).\({k_p}\)is the dimensionless wave number of the variable thickness plate that does not change with x or the plate thickness and can be represented as \({k_p}^{2}=\frac{{{\omega ^2}}}{{{c_p}^{2}}}\); \({c_p}=\sqrt {E/\rho \left( {1 - {\nu ^2}} \right)}\)is the longitudinal wave group velocity in the plate19. \({k_y}\) is the projection of the elastic wave vector in the y direction.

Substitute Eq. (2) into the equation of Eq. (1). Due to that h(x) is a power function about x, If\(x \geqslant 1\), for the same x, the value of x is much larger than the values of \(\frac{{\partial {S^2}(x)}}{{\partial {x^2}}}\), \(\frac{{\partial {S^3}(x)}}{{\partial {x^3}}}\), and \(\frac{{\partial {S^4}(x)}}{{\partial {x^4}}}\), and the higher-order derivative terms of \(S\left( {x,y} \right)\)are discarded to simplify the derivation process. Similarly, discard the higher-order derivative terms of \(A\left( {x,y} \right)\), the equation can be derived by letting the real part be zero.

$${\left| {\nabla S\left( {x,y} \right)} \right|^4}=\frac{{\rho h{\omega ^2}}}{{D{k_p}^{4}}}$$
(3)

By substituting the phase function \(S\left( {x,y} \right)={S_x}\left( x \right)+{k_y}y\) into the Eq. (3), the following can be obtained after arrangement.

$${\left( {\frac{{\partial {S_x}\left( x \right)}}{{\partial x}}} \right)^4}+\frac{{2{k_y}^{2}}}{{{k_p}^{2}}}{\left( {\frac{{\partial {S_x}\left( x \right)}}{{\partial x}}} \right)^2}+\left( {\frac{{{k_y}^{4}}}{{{k_p}^{4}}} - \frac{{\rho h{\omega ^2}}}{{D{k_p}^{4}}}} \right)=0$$
(4)

There are four roots in Eq. (4), which are expressed as

$$\frac{{\partial {S_x}\left( x \right)}}{{\partial x}}= \pm \frac{{{k_y}}}{{{k_p}}}{\left( { - 1 \pm \frac{{{{12}^{\frac{1}{2}}}{k_p}}}{{{k_y}^{2}h}}} \right)^{\frac{1}{2}}}$$
(5)

Expanding Eq. (5) using Newton’s generalized binomial theorem,

$$\begin{gathered} \frac{{\partial {S_x}_{a}\left( x \right)}}{{\partial x}}=\left( \begin{gathered} {12^{{1 \mathord{\left/ {\vphantom {1 4}} \right. \kern-0pt} 4}}}{k_p}^{{ - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}h{\left( x \right)^{ - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} - \frac{1}{2}{12^{ - {1 \mathord{\left/ {\vphantom {1 4}} \right. \kern-0pt} 4}}}{k_y}^{2}{k_p}^{{ - {3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0pt} 2}}}h{\left( x \right)^{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} \hfill \\ - \frac{1}{8}{12^{ - {3 \mathord{\left/ {\vphantom {3 4}} \right. \kern-0pt} 4}}}{k_y}^{4}{k_p}^{{ - {5 \mathord{\left/ {\vphantom {5 2}} \right. \kern-0pt} 2}}}h{\left( x \right)^{{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0pt} 2}}} \hfill \\ \end{gathered} \right) \hfill \\ \frac{{\partial {S_x}_{b}\left( x \right)}}{{\partial x}}= - \left( \begin{gathered} {12^{{1 \mathord{\left/ {\vphantom {1 4}} \right. \kern-0pt} 4}}}{k_p}^{{ - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}h{\left( x \right)^{ - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} - \frac{1}{2}{12^{ - {1 \mathord{\left/ {\vphantom {1 4}} \right. \kern-0pt} 4}}}{k_y}^{2}{k_p}^{{ - {3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0pt} 2}}}h{\left( x \right)^{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} \hfill \\ - \frac{1}{8}{12^{ - {3 \mathord{\left/ {\vphantom {3 4}} \right. \kern-0pt} 4}}}{k_y}^{4}{k_p}^{{ - {5 \mathord{\left/ {\vphantom {5 2}} \right. \kern-0pt} 2}}}h{\left( x \right)^{{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0pt} 2}}} \hfill \\ \end{gathered} \right) \hfill \\ \frac{{\partial {S_x}_{c}\left( x \right)}}{{\partial x}}=\left( \begin{gathered} i{12^{{1 \mathord{\left/ {\vphantom {1 4}} \right. \kern-0pt} 4}}}{k_p}^{{ - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}h{\left( x \right)^{ - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}+\frac{i}{2}{12^{ - {1 \mathord{\left/ {\vphantom {1 4}} \right. \kern-0pt} 4}}}{k_y}^{2}{k_p}^{{ - {3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0pt} 2}}}h{\left( x \right)^{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} \hfill \\ - \frac{i}{8}{12^{ - {3 \mathord{\left/ {\vphantom {3 4}} \right. \kern-0pt} 4}}}{k_y}^{4}{k_p}^{{ - {5 \mathord{\left/ {\vphantom {5 2}} \right. \kern-0pt} 2}}}h{\left( x \right)^{{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0pt} 2}}} \hfill \\ \end{gathered} \right) \hfill \\ \frac{{\partial {S_x}_{d}\left( x \right)}}{{\partial x}}= - \left( \begin{gathered} i{12^{{1 \mathord{\left/ {\vphantom {1 4}} \right. \kern-0pt} 4}}}{k_p}^{{ - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}h{\left( x \right)^{ - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}+\frac{i}{2}{12^{ - {1 \mathord{\left/ {\vphantom {1 4}} \right. \kern-0pt} 4}}}{k_y}^{2}{k_p}^{{ - {3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0pt} 2}}}h{\left( x \right)^{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}} \hfill \\ - \frac{i}{8}{12^{ - {3 \mathord{\left/ {\vphantom {3 4}} \right. \kern-0pt} 4}}}{k_y}^{4}{k_p}^{{ - {5 \mathord{\left/ {\vphantom {5 2}} \right. \kern-0pt} 2}}}h{\left( x \right)^{{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0pt} 2}}} \hfill \\ \end{gathered} \right) \hfill \\ \end{gathered}$$
(6)

Through integral operation on x in Eq. (6), the four roots Sxa, Sxb, Sxc, and Sxd in the expression of the wave-number function \({S_x}\left( x \right)\) can be obtained, which separately represent propagating waves and dissipative wave in the positive and negative x directions, as shown below:

$$\begin{gathered} {S_{xa}}\left( x \right)=\left( \begin{gathered} {12^{{1 \mathord{\left/ {\vphantom {1 4}} \right. \kern-0pt} 4}}}{k_p}^{{ - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}\int {h{{\left( x \right)}^{ - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}dx} - \frac{1}{2}{12^{ - {1 \mathord{\left/ {\vphantom {1 4}} \right. \kern-0pt} 4}}}{k_y}^{2}{k_p}^{{ - {3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0pt} 2}}}\int {h{{\left( x \right)}^{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}dx} \hfill \\ - \frac{1}{8}{12^{ - {3 \mathord{\left/ {\vphantom {3 4}} \right. \kern-0pt} 4}}}{k_y}^{4}{k_p}^{{ - {5 \mathord{\left/ {\vphantom {5 2}} \right. \kern-0pt} 2}}}\int {h{{\left( x \right)}^{{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0pt} 2}}}dx} \hfill \\ \end{gathered} \right) \hfill \\ {S_{xb}}\left( x \right)= - \left( \begin{gathered} {12^{{1 \mathord{\left/ {\vphantom {1 4}} \right. \kern-0pt} 4}}}{k_p}^{{ - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}\int {h{{\left( x \right)}^{ - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}dx} - \frac{1}{2}{12^{ - {1 \mathord{\left/ {\vphantom {1 4}} \right. \kern-0pt} 4}}}{k_y}^{2}{k_p}^{{ - {3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0pt} 2}}}\int {h{{\left( x \right)}^{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}dx} \hfill \\ - \frac{1}{8}{12^{ - {3 \mathord{\left/ {\vphantom {3 4}} \right. \kern-0pt} 4}}}{k_y}^{4}{k_p}^{{ - {5 \mathord{\left/ {\vphantom {5 2}} \right. \kern-0pt} 2}}}\int {h{{\left( x \right)}^{{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0pt} 2}}}dx} \hfill \\ \end{gathered} \right) \hfill \\ {S_{xc}}\left( x \right)=\left( \begin{gathered} i{12^{{1 \mathord{\left/ {\vphantom {1 4}} \right. \kern-0pt} 4}}}{k_p}^{{ - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}\int {h{{\left( x \right)}^{ - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}dx} +\frac{i}{2}{12^{ - {1 \mathord{\left/ {\vphantom {1 4}} \right. \kern-0pt} 4}}}{k_y}^{2}{k_p}^{{ - {3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0pt} 2}}}\int {h{{\left( x \right)}^{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}dx} \hfill \\ - \frac{i}{8}{12^{ - {3 \mathord{\left/ {\vphantom {3 4}} \right. \kern-0pt} 4}}}{k_y}^{4}{k_p}^{{ - {5 \mathord{\left/ {\vphantom {5 2}} \right. \kern-0pt} 2}}}\int {h{{\left( x \right)}^{{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0pt} 2}}}dx} \hfill \\ \end{gathered} \right) \hfill \\ {S_{xd}}\left( x \right)= - \left( \begin{gathered} i{12^{{1 \mathord{\left/ {\vphantom {1 4}} \right. \kern-0pt} 4}}}{k_p}^{{ - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}\int {h{{\left( x \right)}^{ - {1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}dx} +\frac{i}{2}{12^{ - {1 \mathord{\left/ {\vphantom {1 4}} \right. \kern-0pt} 4}}}{k_y}^{2}{k_p}^{{ - {3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0pt} 2}}}\int {h{{\left( x \right)}^{{1 \mathord{\left/ {\vphantom {1 2}} \right. \kern-0pt} 2}}}dx} \hfill \\ - \frac{i}{8}{12^{ - {3 \mathord{\left/ {\vphantom {3 4}} \right. \kern-0pt} 4}}}{k_y}^{4}{k_p}^{{ - {5 \mathord{\left/ {\vphantom {5 2}} \right. \kern-0pt} 2}}}\int {h{{\left( x \right)}^{{3 \mathord{\left/ {\vphantom {3 2}} \right. \kern-0pt} 2}}}dx} \hfill \\ \end{gathered} \right) \hfill \\ \end{gathered}$$
(7)

The far-field waves oscillating in the space can be attained by substituting \({S_x}_{{a,xb}}\left( x \right)\)into the general solution to lateral displacement; \({S_x}_{{c,xd}}\left( x \right)\)correspond to the attenuated near-field waves. For out-of-plane vibration, the near-field term rapidly attenuates to an infinitesimal quantity with the wave propagation, so the term can be ignored. Under the condition, the wave solution to the plate can be expressed as the superposed far-field plane waves. Therefore, Eq. (2) can be represented as

$$w\left( {x,y,t} \right)=\left( {A\left( x \right){e^{i{k_p}{S_{xa}}\left( x \right)}}+B\left( x \right){e^{i{k_p}{S_{xb}}\left( x \right)}}} \right)\left( {{A_y}{e^{ - i{k_y}y}}+{B_y}{e^{i{k_y}y}}} \right){e^{i\omega t}}$$
(8)

In the above equation \(A(x)\),\(B(x)\),\({A_y}\),\({B_y}\) are the vibration displacement amplitude terms. Because the complex stiffness is used to simulate the hysteretic damping effect of the structure, \({k_p}{S_x}\left( x \right)\)and \({k_y}\) in Eq. (8) are also complex numbers. For low-damping structures that\(\eta \ll 1\), there are\({k_y}={k_{y1}}+i{k_{y2}}={k_{y1}}\left( {1 - \frac{{\eta i}}{4}} \right)\) and \({k_p}{S_x}\left( x \right)={k_p}{S_{xa}}\left( x \right)+i{k_p}{S_{xb}}\left( x \right)={k_p}{S_{xa}}\left( x \right) \cdot \left( {1 - \frac{{\eta i}}{4}} \right)\).

The energy density e in the variable thickness plate is the sum of kinetic energy density V and potential energy density T, in which

$$\begin{gathered} T=\frac{1}{2}\rho h{\left( {\frac{{\partial w}}{{\partial t}}} \right)^2} \hfill \\ V=\frac{1}{2}D{\left( {\frac{{{\partial ^2}w}}{{\partial {x^2}}}} \right)^2}+{\left( {\frac{{{\partial ^2}w}}{{\partial {y^2}}}} \right)^2}+2\nu \frac{{{\partial ^2}w}}{{\partial {x^2}}}\left( {\frac{{{\partial ^2}w}}{{\partial {y^2}}}} \right)+2\nu \left( {1 - \nu } \right){\left( {\frac{{{\partial ^2}w}}{{\partial x\partial y}}} \right)^2} \hfill \\ e=T+V \hfill \\ \end{gathered}$$
(9)

For the entire variable thickness plate, the time-averaged energy density and power flow in x and y directions are as follow.

$$\begin{gathered} \overline {e} =\frac{D}{4}\left\{ \begin{gathered} \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\nu \frac{{{\partial ^2}w}}{{\partial {x^2}}}{\left( {\frac{{{\partial ^2}w}}{{\partial {y^2}}}} \right)^*} \hfill \\ +2\nu \left( {1 - \nu } \right)\frac{{{\partial ^2}w}}{{\partial x\partial y}}{\left( {\frac{{{\partial ^2}w}}{{\partial x\partial y}}} \right)^*}+\frac{{\rho h}}{D}\frac{{\partial w}}{{\partial t}}{\left( {\frac{{\partial w}}{{\partial t}}} \right)^*} \hfill \\ \end{gathered} \right\} \hfill \\ \overline {{{q_x}}} ={Q_x}{\left( {\frac{{\partial w}}{{\partial t}}} \right)^*} - {M_{xx}}{\left( {\frac{{{\partial ^2}w}}{{\partial x\partial t}}} \right)^*} - {M_{xy}}{\left( {\frac{{{\partial ^2}w}}{{\partial y\partial t}}} \right)^*} \hfill \\ \overline {{{q_y}}} ={Q_y}{\left( {\frac{{\partial w}}{{\partial t}}} \right)^*} - {M_{yy}}{\left( {\frac{{{\partial ^2}w}}{{\partial x\partial t}}} \right)^*} - {M_{yx}}{\left( {\frac{{{\partial ^2}w}}{{\partial x\partial t}}} \right)^*} \hfill \\ \end{gathered}$$
(10)

Where Q and M separately denote the bending moment and shear force in the plate elements, as expressed below:

$$\begin{gathered} {M_{xx}}= - D\left( {\frac{{{\partial ^2}w}}{{\partial {x^2}}}+\nu \frac{{{\partial ^2}w}}{{\partial {y^2}}}} \right),{M_{yy}}= - D\left( {\frac{{{\partial ^2}w}}{{\partial {y^2}}}+\nu \frac{{{\partial ^2}w}}{{\partial {x^2}}}} \right) \hfill \\ {M_{xy}}={M_{yx}}= - D\left( {1 - \nu } \right)\frac{{{\partial ^2}w}}{{\partial x\partial y}} \hfill \\ {Q_x}= - D\left( {\frac{{{\partial ^3}w}}{{\partial {x^3}}}+\frac{{{\partial ^3}w}}{{\partial x\partial {y^2}}}} \right),{Q_y}= - D\left( {\frac{{{\partial ^3}w}}{{\partial {y^3}}}+\frac{{{\partial ^3}w}}{{\partial {x^2}\partial y}}} \right) \hfill \\ \end{gathered}$$
(11)

Where \(*\) denotes the complex conjugate; \(\nu\)is the Poisson’s ratio. By substituting the far-field solution (8), the time-averaged energy density and power flow of the system can be acquired. let \({k_x}={k_p}{S_x}\left( x \right)\), the energy density and power flow are averaged in unit wavelength space:

$$\begin{gathered} \left\langle {\overline {e} } \right\rangle =\frac{{{k_x}{k_y}}}{{{\pi ^2}}}\int_{{ - \frac{\pi }{{2{k_y}}}}}^{{\frac{\pi }{{2{k_x}}}}} {\int_{{ - \frac{\pi }{{2{k_x}}}}}^{{\frac{\pi }{{2{k_y}}}}} {\left\langle e \right\rangle } } dydx \hfill \\ \left\langle {\overline {q} } \right\rangle =\frac{{{k_x}{k_y}}}{{{\pi ^2}}}\int_{{ - \frac{\pi }{{2{k_y}}}}}^{{\frac{\pi }{{2{k_x}}}}} {\int_{{ - \frac{\pi }{{2{k_x}}}}}^{{\frac{\pi }{{2{k_y}}}}} {\left\langle q \right\rangle } } dydx \hfill \\ \end{gathered}$$
(12)

The time-and-space-averaged energy density and power flow are

$$\begin{gathered} \left\langle {\overline {e} } \right\rangle =\frac{D}{2}{\left( {\left| {{k_{x1}}^{2}} \right|+\left| {{k_{y1}}^{2}} \right|} \right)^2} \hfill \\ \times \left\{ \begin{gathered} {\left| {A\left( x \right)} \right|^2}{\left| {{A_y}} \right|^2}{e^{ - \frac{\eta }{2}\left( {{k_p}{S_{xa}}\left( x \right)+{k_{y1}}y} \right)}}+{\left| {A\left( x \right)} \right|^2}{\left| {{B_y}} \right|^2}{e^{ - \frac{\eta }{2}\left( {{k_p}{S_{xa}}\left( x \right) - {k_{y1}}y} \right)}} \hfill \\ +{\left| {B\left( x \right)} \right|^2}{\left| {{A_y}} \right|^2}{e^{\frac{\eta }{2}\left( {{k_p}{S_{xa}}\left( x \right) - {k_{y1}}y} \right)}}+{\left| {B\left( x \right)} \right|^2}{\left| {{B_y}} \right|^2}{e^{\frac{\eta }{2}\left( {{k_p}{S_{xa}}\left( x \right)+{k_{y1}}y} \right)}} \hfill \\ \end{gathered} \right\} \hfill \\ \end{gathered}$$
(13)
$$\begin{gathered} \left\langle {\overline {{{q_x}}} } \right\rangle =D{k_x}\frac{\omega }{2}\left( {2\left( {{k_{x1}}^{2}+{k_{y1}}^{2}} \right)} \right) \hfill \\ \times \left\{ \begin{gathered} {\left| {A\left( x \right)} \right|^2}{\left| {{A_y}} \right|^2}{e^{ - \frac{\eta }{2}\left( {{k_p}{S_{xa}}\left( x \right)+{k_{y1}}y} \right)}}+{\left| {A\left( x \right)} \right|^2}{\left| {{B_y}} \right|^2}{e^{ - \frac{\eta }{2}\left( {{k_p}{S_{xa}}\left( x \right) - {k_{y1}}y} \right)}} \hfill \\ - {\left| {B\left( x \right)} \right|^2}{\left| {{A_y}} \right|^2}{e^{\eta \left( {{k_p}{S_{xa}}\left( x \right) - {k_{y1}}y} \right)}} - {\left| {B\left( x \right)} \right|^2}{\left| {{B_y}} \right|^2}{e^{\eta \left( {{k_p}{S_{xa}}\left( x \right)+{k_{y1}}y} \right)}} \hfill \\ \end{gathered} \right\} \hfill \\ \left\langle {\overline {{{q_y}}} } \right\rangle =D{k_y}\frac{\omega }{2}\left( {2\left( {{k_{x1}}^{2}+{k_{y1}}^{2}} \right)} \right) \hfill \\ \times \left\{ \begin{gathered} {\left| {A\left( x \right)} \right|^2}{\left| {{A_y}} \right|^2}{e^{ - \frac{\eta }{2}\left( {{k_p}{S_{xa}}\left( x \right)+{k_{y1}}y} \right)}} - {\left| {A\left( x \right)} \right|^2}{\left| {{B_y}} \right|^2}{e^{ - \frac{\eta }{2}\left( {{k_p}{S_{xa}}\left( x \right) - {k_{y1}}y} \right)}} \hfill \\ +{\left| {B\left( x \right)} \right|^2}{\left| {{A_y}} \right|^2}{e^{\frac{\eta }{2}\left( {{k_p}{S_{xa}}\left( x \right) - {k_{y1}}y} \right)}} - {\left| {B\left( x \right)} \right|^2}{\left| {{B_y}} \right|^2}{e^{\frac{\eta }{2}\left( {{k_p}{S_{xa}}\left( x \right)+{k_{y1}}y} \right)}} \hfill \\ \end{gathered} \right\} \hfill \\ \end{gathered}$$
(14)

By comparing Eqs. (13) and (14), the energy intensity and energy density have the following relationship:

$$\left\langle {\overline {q} } \right\rangle = - Ch\left( x \right)\left( {\frac{{\partial \left\langle {\overline {e} } \right\rangle }}{{\partial x}}+\frac{{\partial \left\langle {\overline {e} } \right\rangle }}{{\partial y}}} \right)$$
(15)

Where\(C= - \frac{{4\omega {k_y}^{2}}}{{\eta {{12}^{\frac{1}{2}}}{k_p}}}\)

The Eq. (15) is the constitutive relationship for power-law variable thickness plate. For a conservative system, the energy input, energy dissipation, and energy conduction at a steady state have the following balance relationship. The energy balance is expressed as:

$${\pi _{in}}={\pi _{diss}}+\nabla \left\langle {\overline {q} } \right\rangle$$
(16)

\({\pi _{in}}\) in Eq. (16) indicates the input power density in the system; \({\pi _{diss}}=\eta \omega \left\langle {\bar {e}} \right\rangle\) denotes the dissipation power density of the system. Substitute the Eq. (15) to Eq. (16), the NEFEM constitutive equation for the variable thickness plate can be derived:

$$Ch\left( x \right)\left( {\frac{{\partial ^{2} \left\langle {\overline{e} } \right\rangle }}{{\partial x^{2} }} + \frac{{\partial ^{2} \left\langle {\overline{e} } \right\rangle }}{{\partial y^{2} }}} \right) + Ch^{'} \left( x \right)\left( {\frac{{\partial \left\langle {\overline{e} } \right\rangle }}{{\partial x}} + \frac{{\partial \left\langle {\overline{e} } \right\rangle }}{{\partial y}}} \right) + \eta \omega \left\langle {\overline{e} } \right\rangle = \pi _{{in}}$$
(17)

Where \({\pi _{in}}\)is the external input power density. By solving Eq. (17), the energy density distribution can be obtained. The equation demonstrates the constitutive relationship between input power density and energy density.

Solution of the NEFEM constitutive equation for the variable thickness plate

For Eq. (17), Galerkin weighted residual20 can be used to solve the arithmetic solution to the NEFEM constitutive equation in the sense of weighted integral. The both sides of the constitutive equation are multiplied by the weight function and integral operation is performed in elements, thus attaining the weak variation formulation of the Eq. 

$$\begin{gathered} \int\limits_{\Omega } N ^{T} Ch\left( x \right)\left( {\frac{{\partial ^{2} \left\langle {\overline{e} } \right\rangle }}{{\partial x^{2} }} + \frac{{\partial ^{2} \left\langle {\overline{e} } \right\rangle }}{{\partial y^{2} }}} \right)dxdy + \int\limits_{\Omega } N ^{T} Ch^{'} \left( x \right)\left( {\frac{{\partial \left\langle {\overline{e} } \right\rangle }}{{\partial x}} + \frac{{\partial \left\langle {\overline{e} } \right\rangle }}{{\partial y}}} \right)dxdy \hfill \\ + \int\limits_{\Omega } N ^{T} \eta \omega \left\langle {\overline{e} } \right\rangle dxdy = \int\limits_{\Omega } N ^{T} \pi _{{in}} dxdy \hfill \\ \end{gathered}$$
(18)

By expanding Eq. (18) using Green formula, the following is obtained:

$$\begin{gathered} \int\limits_{\Gamma } N ^{T} \left( {n \times Ch\left( x \right)\nabla \left\langle {\overline{e} } \right\rangle } \right)d\Gamma - \int\limits_{\Omega } {C\nabla \left( {N^{T} \times h\left( x \right)} \right)} \nabla \left\langle {\overline{e} } \right\rangle dxdy \hfill \\ + \int\limits_{\Omega } N ^{T} Ch^{'} \left( x \right)\left( {\frac{{\partial \left\langle {\overline{e} } \right\rangle }}{{\partial x}} + \frac{{\partial \left\langle {\overline{e} } \right\rangle }}{{\partial y}}} \right)dxdy + \int\limits_{\Omega } {\eta \omega N} ^{T} N\left\{ {e^{{\left( e \right)}} } \right\}dxdy - \int\limits_{\Omega } N ^{T} \pi _{{in}} dxdy = 0 \hfill \\ \end{gathered}$$
(19)

After transposition, there is

$$\begin{gathered} {\int\limits_{\Gamma } N ^T}\left( {n \times Ch\left( x \right)\nabla \left\langle {\overline {e} } \right\rangle } \right)d\Gamma - \int\limits_{\Omega } {C \times h\left( x \right)\nabla {N^T}} \nabla \left\langle {\overline {e} } \right\rangle dxdy \hfill \\ +{\int\limits_{\Omega } {\eta \omega N} ^T}\left\langle {\overline {e} } \right\rangle dxdy - {\int\limits_{\Omega } N ^T}{\pi _{in}}dxdy=0 \hfill \\ \end{gathered}$$
(20)

N in Eq. (19) is a weight function, n is the unit normal component \(n=\left( {\cos \theta ,\sin \theta } \right)\) of elements vertical to the boundary \(\Gamma\) and it represents the included angle between the exterior normal of the elements and the x axis; \(\Gamma\)and \(\Omega\)separately denote the element boundary and the entire planar ___domain.

In the elements, the energy density is interpolated through use of the Lagrange interpolation function.

$$\left\langle {\overline {e} } \right\rangle =\sum\limits_{{i=1}}^{n} {{e_i}} {N_i}$$
(21)

\({e_i}\) in Eq. (21) represents the value of energy density at element nodes. By substituting \({e_i}\) in Eq. (20), the following is obtained after simplification:

$$\begin{gathered} {\int\limits_{\Gamma } N ^T}\left( {n \times Ch\left( x \right)\nabla N \cdot {e^{\left( e \right)}}} \right)d\Gamma - \int\limits_{\Omega } {C \times h\left( x \right)\nabla {N^T}} \nabla N \cdot {e^{\left( e \right)}}dxdy \hfill \\ +{\int\limits_{\Omega } {\eta \omega N} ^T}N \times {e^{\left( e \right)}}dxdy - {\int\limits_{\Omega } N ^T}{\pi _{in}}dxdy=0 \hfill \\ \end{gathered}$$
(22)

Equation (22) is written as a matrix equation as

$$\left[ {{K^e}} \right]\left\{ {{e^e}} \right\}=\left\{ {{F^e}} \right\}+\left\{ {{Q^e}} \right\}$$
(23)

Where

$$\begin{gathered} {K^e}=\int\limits_{\Omega } {\left( { - C \times h\left( x \right)\nabla {N^T}\nabla N+\eta \omega {N^T}N} \right)} dxdy \hfill \\ {F^e}={\int\limits_{\Omega } N ^T}{\pi _{in}}dxdy \hfill \\ {Q^e}= - {\int\limits_{\Gamma } N ^T}\left( {n \times Ch\left( x \right)\nabla N \times {e^{\left( e \right)}}} \right)d\Gamma = - {\int\limits_{\Gamma } N ^T}\left( {n \times \left\langle {\overline {q} } \right\rangle } \right)d\Gamma \hfill \\ \end{gathered}$$
(24)

\({K^e}\) is the element stiffness matrix; \({F^e}\) is the column vector of excitation power density input to the element nodes; \({Q^e}\) is the energy flow at the boundary; \(\left\{ {{e^e}} \right\}\) is the energy density. By assembling the matrix equation of all elements and considering the boundary condition, the energy density at each node can be acquired21,22.

Validation and examples

Validation

A variable thickness plate measuring \(\:2\text{m}\times\:2\)m is taken as the research object. The variation of h with x is given by \(h\left( x \right)=0.02*{\left( {1 - {x \mathord{\left/ {\vphantom {x {2.0}}} \right. \kern-0pt} {2.0}}} \right)^2}\). The boundary conditions are shown in the Fig. 2. The power density \(\:1\times\:{10}^{-3}\)W is excited at the center of plate. The material parameters of the variable thickness plate are as follow: Young’s modulus E is \(\:9.6\times\:{10}^{10}\) Pa, Poisson’s ratio υ is 0.35, damping factor η is 0.02, and plate density ρ is \(\:4.62\times\:{10}^{3}\)kg/m3. The plate is meshed. Then 100 four-node variable thickness plate elements (\(\:0.02\text{m}\times\:0.02\text{m}\) for each element) and 121 nodes are generated. The mesh model of the variable thickness plate is illustrated in Fig. 3.

Fig. 2
figure 2

Schematic diagram for the variable thickness plate.

Full size image
Fig. 3
figure 3

NEFEM mesh model of the variable thickness plate.

Full size image

According to the NEFEM established for the variable thickness plate, the distribution of energy density is got. In order to observe the trend of the energy density distribution, energy density levels (EDL) at every nodes of the plate can be calculated according to following equation23.

$$EDL{\text{=}}Lg\left[ {\{ {e^e}\} /1 \times {{10}^{ - 12}}} \right]$$
(25)

In order to validate the NEFEM presented, FEM is applied to predict the dynamic response of a variable thickness plate with the same geometrical parameters, material parameters and boundary conditions (Figs. 4,5). The FEM mesh model is shown in Fig. 6, There are 103,379 nodes, 308,743 tetrahedral elements.

Fig. 4
figure 4

The FEM mesh model for upper surface of the variable thickness plate.

Full size image
Fig. 5
figure 5

The FEM mesh model for lateral surface of the variable thickness plate.

Full size image

Set the excitation frequency as 5000 Hz, the energy density level curve on Y = 1 from NEFEM and FEM are shown in Fig. 6.

Fig. 6
figure 6

The distribution of energy density on Y = 1 from NEFEM and FEM.

Full size image

The trend of energy density level from NEFEM is relatively consistent with that from FEM. It is shown that the NEFEM is accurate.

From the view of calculation time, The NEFEM calculation took 20 s, while the FEM calculation took 229 s. The computational speed of NEFEM is much greater than that of FEM.

It can be seen that the NEFEM can predict accurately dynamic response in less time.

Example

Comparison of results of variable thickness plate with that of uniform plate

To further study if the NEFEM model can show the response characteristic of variable thickness plate, a comparison dynamic response of variable thickness plate with the uniform-thickness plate is analysed. The dynamic response of variable thickness plate calculated by the NEFEM model presented in this paper, and the dynamic response of the uniform-thickness plate is calculated by traditional NEFEM.

The uniform-thickness plate has uniform thickness of 0.02 m, and the length and broad, material parameters, boundary conditions, are same as that of variable thickness plate. The number of elements and nodes are same as that of variable thickness plate. The calculation results are shown in Figs. 7 and 8.

Fig. 7
figure 7

Energy density level distribution of the uniform-thickness plate.

Full size image
Fig. 8
figure 8

Energy density level distribution of the variable thickness plate.

Full size image

As shown in Fig. 8, the energy density level at the excitation point of the variable thickness plate is largest, while energy density level at other nodes decreases with the growing distance to the excitation node. This is because the energy of elastic waves propagates around from the excitation point, which conforms to people’s objective understanding. Because the variable thickness plate is simply supported in two sides y = 0 and y = 2, the energy density level at nodes on the simply supported boundaries is 0. The energy density level is distributed symmetrically about y = 1. Along the direction with variable thickness (x axis), the energy density level at the thin side is higher than that at the thick side. This is because the acoustic black hole effect appears on the plate along the direction of variable thickness. The effect significantly concentrates the vibration energy. As a result, the energy density level is asymmetrically distributed along the x axis.

It can be seen from Fig. 7 that the energy density level is symmetrically distributed in both x and y directions on the uniform-thickness plate. Comparison with Fig. 8 reveals that the NEFEM model can effectively consider influences of the variable thickness when analyzing the high-frequency responses.

Comparison of NEFEM with AEFEM

In order to demonstrate the advantages of the NEFEM model presented in this paper, a comparison NEFEM model with AEFEM model is investigated.

The AEFEM approach the variable thickness plates with large number of constant thickness element. The mesh of AEFEM for the variable thickness plate is shown in Fig. 9. The input power density \(\:1\times\:{10}^{-3}\)W/m2 is applied at the excitation ___location. The excitation frequency is 5000 Hz. And a fixed constraint is applied at the thick side. The rest parameters remain unchanged.

Fig. 9
figure 9

The mesh of AEFEM for the variable thickness plate.

Full size image
Fig. 10
figure 10

The results from the NEFEM, FEM and AEFEM in different element numbers.

Full size image
Table 2 Comparison of the accuracy of results for NEFEM, FEM and AEFEM.
Full size table

It can be seen from Fig. 10; Table 2. that the more elements there are, the closer to the results from NEFEM the results from AEFEM is. If the element number is same, the results from NEFEM model are closer to the results from FEM. This indicates that for the AEFEM, a huge element is needed and the calculation is time consuming. The NEFEM model presented in this paper is more suitable for the dynamic response of variable thickness plates.

Conclusions

The constitutive model based on energy density is crux to predict the high-frequency dynamics response for power-law variable thickness plate in EFEM. At first, the NEFEM constitutive equation was derived for the variable thickness plate according to the energy balance equation and relationship between energy density and strength. Galerkin weighted residual, Lagrange interpolation function, and Green formula were adopted to solve the constitutive equation. Finally, the reasonability of the NEFEM for the variable thickness plate was validated. Main conclusions include as follow.

According to the comparison of the dynamic results from presented NEFEM and the results from FEM, it is shown that the presented NEFEM is reliable and time saving. The presented NEFEM help the engineer to apply the EFEM to predict the high-frequency dynamic response of variable thickness plates.

Compared to FEM, The NEFEM presented in this paper for variable thickness plate does not require excessive meshing and excessive time. However, the NEFEM for variable thickness plate is applicable to the specified power law variable thickness plate that satisfied \(\eta \ll 1\) and power rate variation \(m \geqslant 2\). The other cross-sectional shapes of variable thickness plates are not included.

The energy density level distribution characteristics of the variable thickness plate were analyzed based on examples. The energy density level is symmetrically distributed along the direction of constant thickness in the variable thickness plate; while it is asymmetrically distributed along the direction of variable thickness. The energy density level at the thin side is higher than that at the thick side. This verifies that the energy concentration effect appears invariable thickness structures, which conforms to the acoustic black hole effect.

The geometric nonlinearities and material nonlinearities of variable thickness plates are not considered. It needs to be addressed in the future research to predict the dynamic response of more general engineering structures.