Vibration Analysis of Typical Flat Plate Ship Structures Through a New Method

Abstract

The structural vibration problem in ships and marine structures is one of the core issues in this field. The typical flat plate is a basic structural unit that constitutes the main structures of the hull, such as the cabin, bow and stern, double bottom, superstructure, etc. Typical plate structures and plate frame structures are the basic structures that researchers focus on. This article briefly reviewed the development history of vibration prediction methods for typical flat plate structures and carried out work inspired by different methods and ideas. In order to provide a new tool for predicting the vibration wave energy distribution of typical flat plate structures, a new method for predicting the vibration wave energy distribution of a simple flat plate structure is proposed. This method combines the ideas of ray tracing and the finite element method, discretizes the structure into elements, and represents the transmission of vibration in the structure in the form of energy function mapping between elements and boundaries. After sufficient mapping, the steady-state solution of structural vibration energy density can be obtained. Compared with the finite element method and ray tracing method, this method has better computational efficiency; compared with statistical energy analysis, this method has better spatial resolution.

1. Introduction

In the field of structural mechanics, the problems of the vibration of simple flat plates under a certain load have always been interesting and continuously concerned issues. Simple flat plates are widely used in shipbuilding and ocean engineering, and their own vibration problems are also the “basic problems” that constitute the overall vibration problems of the systems. Generally speaking, the hull structure can be regarded as a complex thin plate structure. The main structures of the hull, such as the cabin, bow and stern, double bottom, superstructure, etc., can be regarded as a composite structure composed of thin plates (elastic plates) connected in different forms. The vibration of the above structures can be refined into an elastic vibration problem of thin plate structures under different boundary conditions. Therefore, the structural vibration characteristics of simple flat plates have always been a concern of designers. On the other hand, simple flat plate problems have been continuously studied by researchers in different ages, and a wealth of data has been accumulated. In the development of new theories, new methods, and new technologies, these data can play the role of verifier and provide a reliable “touchstone” for researchers. It can be seen that vibration problems in simple flat plate structures have high research value. Before introducing the specific work of this article, the author will briefly review the evolution of research methods in the field of structural vibration.

A vast range of numerical methods have been developed for solving noise and vibration problems in mechanical structures. Popular tools include finite element methods, finite volume methods, boundary element methods, and various spectral methods. The core idea of this type of method is to discretize the entire mechanical system or the boundary of the mechanical system into meshes. By establishing and solving a large number of structural dynamics equations, the prediction of unknown response parameters at each node of the system is completed. At low frequencies, the prediction of the characteristics and distribution can be quite satisfactory. However, as the analysis frequency moves to the middle and high frequency bands, the size of the mesh will rapidly decrease with the shortening of the wavelength, followed by mesh number expansion. For large and complex mechanical systems, this means a huge number of calculations and more sensitive solutions to calculation parameters (boundary conditions, damping, mesh quality, etc.) As a result, any slight uncertainty can cause the solution to fail [1].

In order to remove this disadvantage, people have tried to develop a new solution from another perspective—an energy perspective. In 1962 Lyon [2] proposed the power flow theory of bilinearly coupled oscillators and pointed out in 1967 that energy would flow along a gradient of energy density in a mechanical system, just like thermal energy flows along a temperature gradient [3]. This idea has given people a lot of inspiration: if the mechanical system is regarded as a thermodynamic system, the steady-state response solution problem in classical structural mechanics will be transformed into an energy balance problem similar to thermodynamics, which can overcome the above-mentioned problem’s defects, achieve better solution accuracy, and achieve efficiency in the middle- and high-frequency bands. Based on this idea, people have carried out research on energy methods along two main lines [1]: One is based on the assumption of modal density and energy transfer; an energy balance equation including the average input power, internal loss power, and transmission loss power of frequencies in similar modes is established. Fahy [4] in 1987 and Besheara [5] in 1996 proposed the power flow theory of non-conservative strongly coupled oscillators. Together, they formed the basis for later statistical energy analysis (SEA). The other one is based on the wave theory. In the solution domain, the energy balance equation is established, which includes the same waveform in space and frequency as the average input power, loss power, and divergence of the vibration intensity in the domain. After discretization, a linear equation similar to the finite element form can be obtained. In 1989, Nefske and Sung [6] first proposed a power and energy flow mode similar to heat conduction. This method was gradually developed into an energy finite element analysis (EFEA) by combining with finite element technology. Today, the energy method, represented by statistical energy analysis and energy finite element analysis, has become an important tool in the field of mechanical engineering. Among these, statistical energy analysis has been studied more deeply due to its earlier beginnings, and its application depth and breadth in the engineering field also exceeds those of energy finite element analysis.

Statistical energy analysis (SEA) has found widespread application in the automotive and aviation industries, as well as in architectural acoustics. The validity of SEA is based on various assumptions which can be summarized in the following way: (1) It is assumed that subsystems have no memory; that is, coupling constants depend on the properties of subsystems alone. (2) The eigenfunctions of (uncoupled) subsystems are expected to be locally described in terms of random Gaussian fields (‘diffusive wave fields’). These two key assumptions are expected to be valid only in the high-frequency region, for low absorption, and for weakly coupled subsystems having ‘irregular’ shape. The validity of these assumptions is often hard to control, compromising the predictive capability of SEA.

In addition to the energy method, there is another method that is very similar in nature to the energy method but has a large difference in form, which is the ray tracing method (RTM). The idea of the ray tracing method is derived from geometric optics. Energy will propagate through the structure in the form of rays. The energy density at a specific point will be determined by the sum of the ray contributions from all paths from the origin to that point. The entire process follows the law of optics. It is obvious that this method takes into account the flow of all rays. The ray tracing method was originally applied to graphic imaging technology [7]. After people found that the propagation of high-frequency waves and light propagation have great similarities, the ray tracing method began to be widely used in indoor acoustics [8], wireless communication [9], seismic wave prediction and other fields [10]. In the field of mechanical engineering, references [11] described in detail the application of ray tracing to vibration analysis of thin plates. However, tracking rays, including multiple reflections on boundaries, can become cumbersome, in particular when curved surfaces are involved and when including mode conversion at interfaces.

The ray tracing method has many advantages that the statistical energy method does not have. For example, it can better deal with the problem of waves—the number of effective reflections on the interface is relatively small. It can also give more detailed wave energy density distribution and higher spatial resolution, and has better applicability for geometric shapes that are complex and have high degrees of coupling. Comparing the two, the statistical energy method has a higher solution efficiency, but its accuracy depends on the division of the subsystems and the degree of coupling between the subsystems; the ray tracing method has higher solution accuracy—or higher space resolution—and has better adaptability in weird geometry, but its high accuracy depends on having a massive total number of rays, which is detrimental to the efficiency of the solution.

In the past ten years or so, with the development of computer technology, scholars have developed many semi-analytic algorithms and deep learning algorithms to pursue higher calculation accuracy and efficiency. Wang [12,13,14,15] developed a semi-analytical method for the study of functionally graded materials structures and laminated cracked plates; Jin [16,17,18,19] proposed the isogeometric method and achieved good results in vibration research on thin shells, cylindrical shells, and multi-layer structures; Li [20,21,22,23] used the semi-analytic method to study vibration problems in beam–plate–shell structures and the vibro-acoustic response of conical–cylindrical–spherical combined shells, which took into account the effect of free surfaces. Qu [24,25] proposed a general higher-order shear deformation zig-zag theory to analyze composite laminated panels; Tornabene [26,27,28] proposed a higher-order mathematical formulation based on higher-order shear deformation theories, and studied several plate structures such as functional gradient plates and carbon nano plates; M. Arefi [29,30,31] carried out free vibration and buckling analysis on different forms of nano plates based on nonlocal modified higher-order sinusoidal shear deformation theory and nonlocal strain gradient theory. Chong Wang et al. studied the thermal and mechanical properties of composite panels and foam panels [32,33,34]; Chao Yang et al. analyzed the dynamic response of honeycomb composite panels and obtained valuable conclusions [35,36]. Haichao Li et al. used the semi-analytical method to study the dynamic characteristics of plate and shell structures and achieved good results [37,38,39].

The initial motivation of this research comes from the summary and thinking of the above-mentioned methods: Is it possible to develop a new method that not only preserves the essence of the ray tracing method, but also avoids certain assumptions of statistical energy analysis, and can use the existing finite element meshes without introducing too much calculation? Under the guidance of this idea, we began to attempt research on new methods. For a certain mechanical system, by dividing the system into finite elements, that is, discrete elements, and introducing the concept of rays to represent the mechanical energy in the elements, the propagation of energy between the elements is reconstructed into the rays inside the elements and in between the elements. Between reflection and refraction—of course, this follows the laws of optics, and energy is expressed as a function of energy density at the element boundary—through this operation, the calculation weight of the original tens of thousands of rays will be converted into the calculation weight of the energy density function on several element boundaries, thereby greatly reducing the calculation load. In contrast, the smallest calculation element here, that is, the element boundary, is equivalent to the ray in the ray tracing method, or the mesh node in the finite element method, or the subsystem in the statistical energy method. By artificially controlling the mesh size, the accuracy and efficiency of the solution can be flexibly balanced. At the same time, it inherits the advantages of the finite element method and greatly broadens its applicability to complex geometries. The main innovations of this method are shown in Section 1.1.

The structure of the rest of this article is as follows: In the second part, we will introduce the thought and specific process of the combination of ray tracing and the finite element method in detail; after that, we will give the expressions of the energy density of the boundary and the energy density of the receiving point; in the third part, a square plate is taken as an example to predict energy distribution, and the result is compared with the classical solution. Finally, we will briefly summarize the method.

1.1. Key Innovations of the Proposed Method

This method combines the ideas of ray tracing and the finite element method, discretizes the structure into elements (rather than subsystems), and represents the transmission of vibrations in the structure in the form of an energy function (rather than rays) mapped between elements and boundaries. After sufficient mapping, the steady-state solution of the structural vibration energy density can be obtained by integration. The combination of the two is that the discretization of the structure can be directly meshed by the finite element method, and the representation of vibration energy is achieved by the energy function rather than rays. The energy function propagates between the element and the boundary according to the law of optics, and the energy density of each element can be obtained after reaching steady state. This method has the following characteristics:

1. The finite element model of the structure can be directly used for calculation without repeated modeling, thus improving efficiency.

2. High-frequency structural responses are very sensitive to small changes in material properties, boundary conditions, etc. When the wavelength is several orders of magnitude smaller than the typical size of the computational model, the finite element method becomes quite inefficient. The computational efficiency of this method does not change significantly with changes in the analysis frequency.

3. In the ray tracing method, energy propagates in space in the form of rays following the laws of optics. Therefore, it is suitable for calculating the energy distribution of sound fields in complex spaces, and the results have good spatial resolution. However, high-precision results rely on a large number of rays being used in the calculation. Unlike the ray tracing method, this method represents the energy function with orthogonal polynomials, and its minimum calculation unit is the element boundary mapping, which is much smaller than the number of rays, the minimum calculation unit in the ray tracing method.

4. The computational accuracy of statistical energy analysis depends on appropriate subsystem partitioning, and there is no spatial resolution for the energy distribution within the subsystem. A major advantage of this method is that the choice of subsystem partitioning is no longer critical, because the SEA assumption related to this has been removed. The meshing of this method is more convenient, and better spatial resolution can be obtained as needed.

5. This method can be regarded as a complex form of statistical energy analysis in terms of geometry, and as a simple form of ray tracing in terms of algorithms. The first two methods can be roughly regarded as two limit states of this method. According to actual needs, the number of elements and the parameters of the energy function can be manually adjusted to balance calculation accuracy and efficiency. The finite element form of mesh division provides this method with better applicability to complex structures.

2. Thought Process and Theoretical Derivation

The method introduced in this article is derived from the ray tracing method and finite element discrete technology. It uses the density function to represent the energy distribution of the ray on the boundary. Each boundary in space is the minimum calculation element. Multiple energy mappings between simulations simulate multiple reflections or refractions of rays between different boundaries in space. By dividing a complex system into several elements, the energy density on the boundary of each element is mapped to other boundaries according to a certain mapping relationship, so as to simulate the reflection on the element boundary or the propagation of refracted rays on the boundaries of adjacent elements.

First, the system is discretized into several elements, and the element where the source point is located is determined. According to a certain mapping relationship, the energy of the source point is mapped to the boundary of the element to obtain the initial boundary energy density. The initial boundary energy density of each edge will be mapped to the other boundaries of the element. At the same time, the initial energy density on other boundaries will also be mapped. The actual structure usually consists of many elements, so the energy densities on the element boundaries will not only be mapped to each other within the element, but also to the boundaries of adjacent elements. When the energy of the ray is reduced to a negligible level after going through several mappings, the density obtained by all mappings is superimposed to obtain the energy density after multiple reflections on each side. Finally, according to the position of the receiving point in the actual structure, the element is determined based on location, and the energy density functions of the boundaries of the element are mapped to the location of the receiving point to obtain the energy density of the receiving point.

2.1. Initial Mapping

Initial mapping is the first step, and its main purpose is to convert the external excitation load input into the initial boundary energy density. Generally speaking, the loading method for concentrated excitations is used to simplify it to a point source and act on the centroid of the selected loading element (source element). As shown in Figure 1, the source element Ω1 is a triangular mesh with a boundary Γ1, which can be discretized into three edges, E11, E12, and E13. It should be noted that E11, E12, and E13 here are the “inner edges” (blue lines) of Ω1, and their normal directions point to the inside of the element. Corresponding to the “inner edges” are the “outer edges” (green lines), the meaning of which will be introduced in the next part. There is a point source inside Ω1 (red point), and the point source radiates energy around. The energy input from the point source will be propagated to E11, E12, and E13 in the form of mapping. We introduce the initial boundary energy density ρ(Γ1) to describe this mapping process. ρ(Γ1) can be disassembled into the sum of ρ(E11), ρ(E12), and ρ(E13).

Here, we introduce the phase space coordinate system $x=(s,p)$, where $s$ represents position and $p$ represents momentum. Let the point source radiate energy at a frequency of ${\omega }_0$ and an intensity of $R_0$; then, the spatial energy at the distance $r$ from the source $r_0$ is:

$\rho \left(r,p\right)={\displaystyle \frac{R_0}{2\pi v_g}}{\displaystyle \frac{e^{-\mu \left|r-r_0\right|}}{\left|r-r_0\right|}}$

Here, $\mu$ is damping coefficient, and $v_g$ is the group velocity. According to the geometric relationship shown in Figure 2, the initial boundary energy density ${\rho }_0$ can be obtained as:

(1)${\rho }_0\left(s,p\right)={\displaystyle \frac{R_0}{2\pi {\omega }_0}}{\displaystyle \frac{e^{-\mu D\left(s,r\right)}}{D\left(s,r\right)}}\delta \left(p-k_0\mathit{sin}\alpha \left(s,r\right)\right)\mathit{cos}\alpha \left(s,r\right)$

Here, $D\left(s,r\right)$ is the distance from the boundary points and the source point. We can now express the initial boundary energy density ${\rho }_0$ in the form of polynomials and a series of unknown coefficients:

(2)${\rho }_0\left(x_0\right)={\displaystyle \sum _n}{f}_{0,n}{F}_n\left(x_0\right),\hspace{1em}\hspace{1em}{f}_{0,n}={\displaystyle \frac{\left({\rho }_0,F_n\right)}{\left(F_n,F_n\right)}}$

Here:

(3)$f_{0,\left(b,\beta \right)}={\displaystyle \frac{R_0}{2\pi {\omega }_0}}\sqrt{{\displaystyle \frac{1}{A_b{k}_b}}}{\displaystyle \frac{2\beta +1}{2}}{\int }_0^L\lambda \left(\alpha \left(s\right)\right)e^{-\mu D\left(s\right)}{P}_{\beta }\left(\mathit{sin}\alpha \left(s\right)\right){\displaystyle \frac{\mathit{cos}\alpha \left(s\right)}{D\left(s\right)}}ds$

It can be obtained that the integral term is only a function of $s$, and the unknown coefficient vector $f_{0,n}$ can be obtained using the numerical integration method.

(4)$f_{0,n}=\left[\begin{array}{cccccc}{f}_{0,\left(\mathrm{1,1}\right)}& f_{0,\left(\mathrm{1,2}\right)}& \cdots & f_{0,\left(1,N\right)}& \begin{array}{cc}{f}_{0,\left(\mathrm{2,1}\right)}& \cdots \end{array}& f_{0,\left(3,N\right)}\end{array}\right]$

2.2. Boundary Mapping

The element boundary is the smallest calculation element in this method, and energy mapping between the boundaries is the most critical step of this method. In this section, a two-dimensional triangular element is taken as an example to explain the whole process of boundary energy mapping in detail. According to the principle of the discrete element in the finite element method, complex structures are usually discretized into multiple triangular elements. Energy is transferred to adjacent elements through the boundary. Among them, adjacent elements usually share a boundary, and the boundary is divided into several half-boundaries, and each half-boundary corresponds to an element. The boundary contains all the rays that have reflected or refracted from this boundary. For a straight boundary, it includes the ray angle α ∈ [−π/2, π/2], and the ray energy transmitted to the inside of the element or other boundaries of the element in any direction. Among them, the ray angle α represents the angle between the ray and the boundary normal. When a boundary is adjacent to two elements, the boundary is divided into two half boundaries, as shown in Figure 3. For element Ω1, edge E1 is the inner boundary and E2 is the outer boundary; for element Ω2, edge E2 is the inner boundary and E1 is the outer boundary.

Here, we explain the considerations involved in using triangular meshes. The use of triangular meshes is essentially carried out to more conveniently simulate the reflection of energy functions between mesh boundaries. When quadrilateral meshes are used, boundary mapping will occur between two non-adjacent sides. Compared with situations where three sides of a triangular mesh are adjacent to each other, quadrilateral meshes will be more complicated in integral calculation. Therefore, this method uses triangular meshes. Considering that all quadrilateral meshes can be divided into two triangular meshes, there is no essential difference between these two division methods, but triangular meshes will be simpler in early preprocessing and calculation.

The coordinate mapping relationship between boundaries is explained in detail below. Figure 4 shows two adjacent triangular elements, where the positive direction of all boundary line coordinates is counterclockwise, as shown by the pointing arrow in Figure 5. A ray on the boundary 1 in element 1 that starts at (s₀, p₀) propagates on boundary 2 of the element at an angle α₀; the incident angle is α₁, and intersects at (s₁, p₁). The rays then form reflections and refractions at (s₁, p₁). At this time, boundary 2 is adjacent to the two elements, so it is divided into two half boundaries, one of which contains the energy reflected back to element 1, while the other half of the boundary contains the energy refracted into element 2.

The actual structure usually contains a large number of elements, and the calculation of the boundary map is huge. Therefore, it is necessary to establish a mapping of the overall phase space, that is, including the mutual mapping between all boundaries. Let the two-dimensional phase space coordinate of the boundary point be X = (s, p), which contains position and momentum information. In a triangle element, let two adjacent half edges be denoted as E₁ and E₂, respectively. The energy density ρ₁(X₁) on the boundary E is mapped onto the boundary E₂ to obtain the energy density ρ₂(X₂). For the case of no energy loss or reflection refraction, ${\rho }_2\left(\Phi \left(X_1\right)\right)={\rho }_1\left(X_1\right)$. Here, the Perron–Frobenius operator is introduced [40]:

(5)$\rho \left(X,\tau \right)={\mathcal{L}}^{\tau }[\rho ](X)=\int \delta \left(X-{\Phi }^{\tau }\left(Y\right)\right)\rho \left(Y,0\right)dY$

The energy map can be expressed as an integral form:

(6)${\rho }_2\left(X_2\right)=\left\{T\rho \right\}\left(X_2\right)=\int \lambda \left(X_2\right)e^{-\mu D\left(X_2,X_1\right)}\delta \left(X_2-\Phi \left(X_1\right)\right)\rho \left(X_1\right)d{X}_1$

Here, $\lambda \left(X_2\right)$ represents the energy loss of the ray in the process of boundary reflection or refraction, which is similar to the concept of the coupling loss factor in the statistical energy method. $e^{-\mu D\left(X_2,X_1\right)}$ Represents energy loss due to the damping factor of the medium during the ray propagation, where μ represents the damping factor, and $D\left(X_2,X_1\right)$ represents the propagation distance from boundary 1 to boundary 2.

At this time, the polynomial and a series of unknown coefficients are used to represent the initial energy density function and the energy density function after mapping:

(7)${\rho }_1\left(X_1\right)=\sum _n{f}_n{F}_n\left(X_1\right),{\rho }_2\left(X_2\right)=\sum _m{g}_m{F}_m\left(X_2\right)$

The basis function of the position space is a piecewise constant function, and the basis function of the momentum space is a Legendre polynomial; that is:

(8)$F_{\left(b,\beta \right)}\left(s,p\right)={\displaystyle \frac{1}{\sqrt{L_b{k}_b}}}{\delta }_b{P}_{\beta }\left({\displaystyle \frac{p}{k_b}}\right)$

where b represents the sequence number of the boundary; β represents the order of the polynomial, $b=\mathrm{1,2},3;\beta =\mathrm{0,1},2,...,N;$L_b represents the length of the boundary b; and ${\delta }_b$ represents that when s ∈ b, ${\delta }_b$ = 1, otherwise 0. Based on the inner product of the polynomial, there will be:

(9)$\left(F_{\left(b_1,{\beta }_1\right)}\left(s_1,p_1\right),F_{\left(b_2,{\beta }_2\right)}\left(s_2,p_2\right)\right)={\displaystyle \frac{2}{2\beta +1}}\delta \left(b_2-b_1\right)\delta \left({\beta }_2-{\beta }_1\right)$

Let $n=\left(b,\beta \right)$, substituting Equation (8) into Equation (7), we obtain:

(10)$\begin{array}{c}{T}_{mn}={\displaystyle \frac{2{\beta }_2+1}{2\sqrt{L_{b_1}{L}_{b_2}{k}_{b_1}{k}_{b_2}}}}\\ ={\int }_0^{L_{b_1}}{\int }_{-k_{b_1}}^{k_{b_1}}{\lambda \left(\Phi \left(X_1\right)\right)e}^{-\mu D\left(s_2,s_1\right)}{\delta }_{b_1}(s_1){\delta }_{b_2}(s_2)P_{{\beta }_1}\left({\displaystyle \frac{p_1}{k_{b_1}}}\right)P_{{\beta }_2}\left({\displaystyle \frac{p_2}{k_{b_2}}}\right)dpds\hfill \end{array}$

The integral relationship is shown in the Figure 6. At this time, the integral variables in the operator T are s and p. With the ray angles α and s as variables, the operator can be expressed as:

(11)$\begin{array}{c}{T}_{\mathit{m}\mathit{n}}={\displaystyle \frac{2{\beta }_2+1}{2\sqrt{L_{b_1}{L}_{b_2}}}}\sqrt{{\displaystyle \frac{k_{b_1}}{k_{b_2}}}}\\ {\int }_0^{L_{b_1}}{\int }_{-\pi /2}^{\pi /2}\lambda \left(\Phi \left(X_1\right)\right)e^{-\mu D\left(s_2,s_1\right)}{\delta }_{b_1}\left(s_1\right){\delta }_{b_2}\left(s_2\right)P_{{\beta }_1}\left(\mathit{sin}{\alpha }_2\right)P_{{\beta }_2}\left(\mathit{sin}{\alpha }_1\right)\cos {\alpha }_{1}d{\alpha }_{1}ds\hfill \end{array}$

Having determined the transfer operator T, the stationary phase space density ${\rho }_{\infty }$ on the boundaries of each of the mesh cells is obtained by iteration; that is [41]:

(12)${\rho }_{\infty }=\sum _{n=0}^{\infty }{T}^n{\rho }_0={(1-T)}^{-1}{\rho }_0$

Here, ${\rho }_0$ is the initial source density projected onto the mesh boundaries that will be discussed in Section 2.2. $n$ represents the number of times boundary mapping occurs, and when $n=0$, no boundary mapping occurs, which means that energy only exists in the internal boundary of the cell where the source point is located. When $n=1$, one time boundary mapping occurs. The energy on one boundary is only transferred to the two adjacent half-boundaries, which correspond to the two adjacent edges of the triangle, as shown in Figure 3. Therefore, one time mapping can only be transferred to three half-boundaries, as shown as Figure 2.

As shown in Figure 2, when the boundary mapping does not occur, $n=0$, and the energy from the source point is only mapped to the three half-boundaries (red lines) inside the cell. After one time boundary mapping, $n=1$, the energy is mapped to the adjacent three half-boundaries (yellow lines). If the mapping process of $n=1$ is disassembled, it can be regarded as mapping $E_{11}$ to $E_{21}$ and $E_{31}$; then, $E_{11}$ corresponds to phase space $X$, and $E_{21}+E_{31}$ corresponds to phase space $X^{\prime }$; that is, the energy on $E_{11}$ should be equal to the sum of the energy on $E_{21}$ and $E_{31}$, which means the following relationship exists:

(13)$\int {\varrho }_{E_{11}}\left(X\right)\mathrm{d}X=\int ({\varrho }_{E_{21}-E_{11}}\left(X^{\prime }\right)+{\varrho }_{E_{31}-E_{11}}\left(X^{\prime }\right))\mathrm{d}{X}^{\prime }$

The same relationship applies to $E_{12}$ and $E_{13}$:

(14)$\int {\varrho }_{E_{12}}\left(X\right)\mathrm{d}X=\int ({\varrho }_{E_{31}-E_{12}}\left(X^{\prime }\right)+{\varrho }_{E_{41}-E_{12}}\left(X^{\prime }\right))\mathrm{d}{X}^{\prime }$

(15)$\int {\varrho }_{E_{13}}\left(X\right)\mathrm{d}X=\int ({\varrho }_{E_{41}-E_{13}}\left(X^{\prime }\right)+{\varrho }_{E_{21}-E_{13}}\left(X^{\prime }\right))\mathrm{d}{X}^{\prime }$

Combine Formulas (13)–(15) to obtain:

(16)$\int \left({\varrho }_{E_{11}}\left(X\right)+{\varrho }_{E_{12}}\left(X\right)+{\varrho }_{E_{13}}\left(X\right)\right)\mathrm{d}X=\int \left({\varrho }_{E_{21}}\left(X^{\prime }\right)+{\varrho }_{E_{31}}\left(X^{\prime }\right)+{\varrho }_{E_{41}}\left(X^{\prime }\right)\right)\mathrm{d}{X}^{\prime }$

In Formula (16), $E_{11}+E_{12}+E_{13}$ corresponds to the phase space $X$, and $E_{21}+E_{31}+E_{41}$ corresponds to the phase space $X^{\prime }$. In this way, substituting $n=1$ into Formula (12), we can obtain:

(17)${\varrho }_1=\sum _{n=0}^1{T}^n{\varrho }_0=I{\varrho }_0+T^1{\varrho }_0$

Here, we obtain:

(18)$I{\varrho }_0={\varrho }_{E_{11}}+{\varrho }_{E_{12}}+{\varrho }_{E_{13}}$

(19)$T^1{\varrho }_0={\varrho }_{E_{21}}+{\varrho }_{E_{31}}+{\varrho }_{E_{41}}$

Therefore, after one time boundary mapping, the following formula holds:

(20)${\varrho }_1={\varrho }_{E_{11}}+{\varrho }_{E_{12}}+{\varrho }_{E_{13}}+{\varrho }_{E_{21}}+{\varrho }_{E_{31}}+{\varrho }_{E_{41}}$

If we look at Formulas (19) and (22) together, the following problem arises:

(21)${\varrho }_{E_{21}}+{\varrho }_{E_{31}}+{\varrho }_{E_{41}}=T^1{\varrho }_0$

(22)${\varrho }^{\prime }\left(X^{\prime }\right)=\left\{T\varrho \right\}\left(X^{\prime }\right)$

As shown in Figure 7, boundary mapping is the core part, and it is also the part that best reflects the core of the method. In fact, we can make this understanding intuitively: $x=(s,p)$ represents the phase space coordinates, composed of all the coordinates $s$ of the boundary and all the directions $p$ of the ray emitted by each coordinate; $\varrho \left(x\right)$ represents the boundary energy density of the mapping edge; and ${\varrho }^{\prime }\left(x^{\prime }\right)$ represents the boundary energy density of the receiving edge. For ease of understanding, one may assume that there is a light source that can slide along the boundary with the angle range of the light $[-{\displaystyle \frac{\pi }{2}},{\displaystyle \frac{\pi }{2}}]$. In this way, all possible positions of the light source and the direction of the emitted light constitute all $x=(s,p)$. If a fixed $x=(s,p)$, then $\varrho \left(x\right)$ represents the intensity of light emitted by this light source at this position and direction.

In addition, regarding the boundary mapping part, we will give a simple example in Appendix A to further illustrate it.

2.3. Final Mapping

After completing the initial mapping and boundary mapping, all the boundaries of all elements have corresponding steady-state energy densities. The final mapping procedure is the remapping of the inner boundary energy density belonging to each element back to the inside of the element, so as to obtain the state energy density of each element. This process can be seen in Figure 8.

The sound source point or energy source point is mapped through the initial boundary to obtain the equivalent initial energy density on the element boundary, and then the Perron–Frobenius operator is mapped to obtain the energy density function on each element boundary. Finally, by mapping the energy density function on the boundary to the inside of the element, the energy density at any point in the element can be obtained.

(23)$\begin{array}{c}{\rho }_{sp}\left(r\right)={\int }_0^{\infty }{\int }_0^{2\pi }{e}^{-\mu \left|r-r_{bd}\right|}\rho \left(s\left(r,\phi \right),k\sin \alpha \left(r,\phi \right)\right){\displaystyle \frac{{\omega }_0}{v_{g,0}}}\delta \left(k-k_0\right)d\phi dk\\ ={\displaystyle \frac{{\omega }_0{k}_0}{v_{g,0}}}{\int }_0^{2\pi }{e}^{-\mu \left|r-r_{bd}\right|}\rho \left(s\left(r,\phi \right),k\sin \alpha \left(r,\phi \right)\right)d\phi \hfill \end{array}$

where φ is the angle describing the direction of ray propagation in the polar coordinate system. $s\left(r,\phi \right)$ and $\alpha \left(r,\phi \right)$ represent the position of the intersection of the ray and the boundary and the angle between the ray and the boundary normal, respectively. The integral variable in Equation (23) is changed to the boundary position coordinate s, which can be obtained as follows:

(24)${\rho }_{\mathrm{sp}}\left(r\right)={\displaystyle \frac{{\omega }_0{k}_0}{v_{g,0}}}{\int }_0^{2\pi }{e}^{-\mu D\left(s,r\right)}\rho \left(s,k_0\mathit{sin}\alpha \left(s,r\right)\right){\displaystyle \frac{\mathit{cos}\alpha \left(s,r\right)}{D\left(s,r\right)}}ds$

where $D\left(s,r\right)$ is the distance between the intersection point on the boundary and the pointr. Extend the base function in the following form:

(25)$\rho \left(s,p\right)=\sum _{b,\beta }{\displaystyle \frac{f_{b,\beta }}{\sqrt{L_b{k}_b}}}{\delta }_b{P}_{\beta }\left({\displaystyle \frac{p}{k_b}}\right)$

Substituting into Formula (24), the space energy density can be written as:

(26)${\rho }_{\mathrm{sp}}\left(r\right)={\displaystyle \frac{{\omega }_0}{v_{g,0}}}\sum _{b,\beta }{f}_{b,\beta }\sqrt{{\displaystyle \frac{k_b}{L_b}}}{\int }_0^{L_b}{e}^{-\mu D\left(s,\mathit{r}\right)}{\delta }_b{P}_{\beta }\left(\mathit{sin}\alpha \left(s,r\right)\right){\displaystyle \frac{\mathit{cos}\alpha \left(s,r\right)}{D\left(s,r\right)}}ds$

Here, we explain the meaning of Formula (26); it can be roughly explained as follows: by mapping the energy density on the internal boundary of the triangular element to the receiving point to obtain the energy density of the point, in this method, each element can obtain an energy density value. This is somewhat similar to the ray tracing method in which rays are gathered to the receiving point. In this method, the energy is discretized into energy density on each element boundary rather than a complete ray, so these distance and angle parameters cannot be directly analogous to the distances and angles in the ray tracing method.

The sum in the formula represents all the boundaries of all the elements containing the point r; that is, the energy of all the boundaries of the elements is mapped to the sum of the points.

2.4. Brief Process

Here is a brief operation process for calculation using this method, shown as Figure 9:

• (a). Select a structure (a simple flat plate, in this article).

• (b). Meshing (triangle elements).

• (c). Select the loading element (the element where the load is located) and use Equations (3) and (4) in Section 2.1 to obtain the boundary energy density of the loading element, which is the initial boundary energy density.

• (d). The initial boundary energy density of each edge will be mapped to the other boundary of the element, and at the same time, it will be mapped to the boundary of the adjacent element. The steady-state value of the boundary energy density can be obtained using Equation (12). Section 2.1 mainly describes this process, in which Equation (10) and Equation (11) can be used to obtain the mapping operator between different boundaries.

• (e). The energy density of each boundary of the element is mapped to the receiving point, and the energy density value of the receiving point can be obtained by Equation (26) in Section 2.3. After obtaining the energy values of all elements, the energy distribution of the entire structure is also obtained.

Figure 9

Brief operation process.

[Figure omitted. See PDF]

3. Single Flat Plate Simulation

After completing the theoretical work, we use the examples in the literature [42] as simulation objects to explore the usability of our method. As shown in Figure 10: The calculation model is a square plate of 1 m × 1 m with simple supports around it. The excitation load acts on the center of the plate with a value of F = 0.01 N, the thickness is 0.001 m, the Young’s modulus E = 7.1 × 10¹⁰ N/m², and the density ρ = 2700 kg/m³.

In this article, four simulation conditions are calculated for the flat plate model, and the “exact solution” given in the literature [42] and the predicted value obtained from the energy finite element analysis (EFEA) are set as comparisons to allow for discussion of the effects of the energy density accuracy (see Table 1 for details).

In the first calculation example, the calculated frequency and damping are f = 239 Hz and η = 0.05; the “exact” solution of the energy distribution along the diagonal of the plate (Exact), the energy finite element analysis prediction value (EFEA), and the predicted value by the method proposed in this article (Numerical) are shown in Figure 11. The energy distribution and flow obtained by the method in this article are also given in Figure 12 and Figure 13 together. It can be seen that the shapes of the three curves are almost completely different. The “exact” solution exhibits rich “gully” characteristics. The EFEA predicted value is the smoothest overall. Compared with the “exact” solution, the predicted value of the method in this article is smoother, and it has more fluctuations compared to EFEA. Generally speaking, the predicted value of EFEA is closer to the “exact” solution than this method, and the predicted value of the method in this article near the point of excitation is closer to the “exact” solution. From the energy distribution map (Figure 12) and energy flow map (Figure 13), the flat plate energy distribution obtained by the method presented in this article presents a certain “hexagonal” feature. This is due to the triangulation method used in this method to mesh the structure, and the energy mapping relationship based on the triangular mesh. This phenomenon can be reduced by increasing the number of elements, but it cannot be eliminated fundamentally.

In the first calculation example, the calculated frequency and damping are f = 239 Hz, η = 0.05; the “exact” solution of the energy distribution along the diagonal of the plate (Exact), the energy finite element analysis prediction value (EFEA), and the predicted value by the method proposed in this article (Numerical) are shown in Figure 14. The energy distribution and flow obtained by the method in this article are also given in Figure 15 and Figure 16 together. It can be seen that as the damping increases, the “slopes” of the three curves increase significantly, which represents the acceleration of the energy attenuation, and that the three curves are smoother than when the damping is low. In general, the method in this article results in values closer to the EFEA prediction value, and it is still closer to the “exact solution” near the excitation point. According to the energy density distribution map (Figure 15) and flow map (Figure 16) of the method in this article, the increase in damping makes both the energy distribution map and the energy flow map smoother, and the “hexagonal” characteristics of the energy flow map are weakened to some extent.

In the third and fourth examples, the calculation frequency is increased to f = 487 Hz, and the damping is still η = 0.05 and η = 0.20. The “exact” solution of the energy distribution along the diagonal of the plate (Exact), the energy finite element analysis prediction value (EFEA), and the predicted value by the method proposed in this article (Numerical) are shown in Figure 17 The energy distribution and flow obtained by the method in this article are also given in Figure 18, Figure 19, Figure 20, Figure 21 and Figure 22, together. It can be seen that with an increase in frequency and damping, the predicted value of EFEA and the predicted value of this method become closer to the “exact” solution, and the three curves also become shinier. On the whole, the prediction value of the method in this article will be closer to the EFEA prediction value under high frequency and high damping conditions. It is lower than EFEA in remote areas far from the excitation point and higher than EFEA near the excitation point.

4. Conclusions

The typical flat plate is the basic structural unit that constitutes the main structures of the hull, such as the cabin, bow and stern, double bottom, superstructure, etc. Therefore, the typical plate vibration problem has attracted extensive attention from researchers. By combining the ideas behind the finite element method and the ray tracing method, this article proposes a new method for solving the vibration energy distribution of a simple flat structure. After meshing the flat structure, the energy flow in the mesh is simulated by using rays to construct an energy map in a single mesh, and it is generalized to adjacent meshes. Finally, the energy density and energy flow direction of each mesh can be solved. This method combines the advantages of the finite element method and the ray tracing method, and can be used as a new attempt to solve the energy distribution in a simple plate. Through previous discussions and research on the theory and with calculation examples, we can temporarily summarize some of the characteristics of the method proposed in this article.

• (1). The basic idea of the method proposed in this article is to divide the structure based on the triangular mesh, comparing its energy flow in the form of energy mapping. Therefore, the specific division of the element will affect the final calculation result by changing the energy mapping function. The aforementioned “hexagonal” feature is a concrete manifestation of this effect, and this effect cannot be completely eliminated.

• (2). The method proposed in this article is still a semi-analytical method. The prediction result is a smooth approximation of the “exact” solution, as shown above, so it cannot reflect the details of the energy distribution in each area of the plate. This feature is similar to EFEA, but compared with EFEA, the method in this article can “slightly” reflect some details.

• (3). With increases in frequency and damping, the method proposed in this article becomes closer to the “exact solution”. In terms of energy distribution, compared with EFEA, in the edge area far from the excitation point, the predicted value of EFEA is higher than the method in this article and lower than the method in the vicinity of the excitation point. Only in the vicinity of the excitation point can the method in this article give prediction results closer to the “exact” solution.

• (4). Based on the above points, the method proposed in this article provides a new idea for the energy prediction of simple flat plates. By further deepening the theoretical algorithm, better prediction results could be obtained.

Here, we can summarize some characteristics of this method: This method borrows the idea of ray tracing in physics, but the difference is that we use energy functions instead of rays to represent the transmission of vibration in the structure; mathematically, it borrows the idea of the finite element method, discretizes the energy function by dividing the structure into elements, and establishes an energy function mapping matrix for the entire structure to simulate the transmission of vibrations in the structure. The calculation accuracy of this method does not depend on a large number of rays and refined modeling. It allows for the adjusting of the parameters of the energy function and the number of elements according to actual needs to balance calculation accuracy and efficiency. It also shows potential in application to complex structures, and can be used as a new way to analyze structural vibration.

Author Contributions

Methodology, J.L. and H.L.; formal analysis, Y.T.; writing—original draft, S.L.; writing—review and editing, Y.T.; supervision, H.L. and X.M.; project administration, X.M. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Data is contained within the article.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

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

Figures and Table

Generate image description

Figure 1. Initial mapping.

Generate image description

Figure 2. Initial mapping of sound source to the boundary in the element. The energy of the source point will be mapped to each boundary to obtain the initial boundary energy density.

Generate image description

Figure 3. Half-boundary diagram. Adjacent elements usually share a boundary, and the boundary is divided into 2 half-boundaries, then each half-boundary corresponds to an element.

Generate image description

Figure 4. Boundary mapping. (a) Boundary mapping for an element. (b) Boundary mapping for an edge.

Generate image description

Figure 5. Ray mapping between boundaries. There are two half-boundaries adjacent to element 1 and element 2. One half-boundary contains the energy reflected back to element 1, and the other half-boundary contains the energy refracted into element 2.

Generate image description

Figure 6. Conversion relationship between angle variable and distance variable. Assuming that [Forumla omitted. See PDF.] is only a function of the angle α, the double integration can be simplified to a single integration with the angle α as the integration variable.

Generate image description

Figure 7. Understandings of boundary mapping.

Generate image description

Figure 8. Final mapping.

Generate image description

Figure 10. Simple support plate. The simulation model is a simple flat structure, surrounded by simple support constraints, and the excitation load acts on the center.

Generate image description

Figure 11. The energy distribution along the diagonal of the plate. Here, the frequency f = 239 Hz and the damping η = 0.05. Exact (blue line) represents the solution given by the literature [42], EFEA (black line) represents the predicted value obtained from the energy finite element analysis, and Numerical (red line) represents the predicted value of the method proposed in this article. The reference value of energy density is 1 × 10−12 J/m2.

Generate image description

Figure 12. The energy distribution predicted by the method proposed in this article. Here, the frequency f = 239 Hz and the damping η = 0.05. The reference value of energy density is 1 × 10−12 J/m2.

Generate image description

Figure 13. The energy flow map of the energy distribution predicted by the method proposed in this article. Here, the frequency f = 239 Hz, and the damping η = 0.05. The reference value of energy flow is 1 × 10−12 W/m2.

Generate image description

Figure 14. The energy distribution along the diagonal of the plate. Here, the frequency f = 239 Hz and the damping η = 0.20. Exact (blue line) represents the solution given by the literature [42], EFEA (black line) represents the predicted value obtained from the energy finite element analysis, and Numerical (red line) represents the predicted value of the method proposed in this article. The reference value of energy density is 1 × 10−12 J/m2.

Generate image description

Figure 15. The energy distribution predicted by the method proposed in this article. Here, the frequency f = 239 Hz, and the damping η = 0.20. The reference value of energy density is 1 × 10−12 J/m2.

Generate image description

Figure 16. The energy flow map of the energy distribution predicted by the method proposed in this article. Here, the frequency f = 239 Hz, the damping η = 0.20. The reference value of energy flow is 1 × 10−12 W/m2.

Generate image description

Figure 17. The energy distribution along the diagonal of the plate. Here, the frequency f = 487 Hz and the damping η = 0.05. Exact (blue line) represents the solution given by the literature [42], EFEA (black line) represents the predicted value obtained from the energy finite element analysis, and Numerical (red line) represents the predicted value of the method proposed in this article. The reference value of energy density is 1 × 10−12 J/m2.

Generate image description

Figure 18. The energy distribution predicted by the method proposed in this article. Here, the frequency f = 487 Hz, and the damping η = 0.05. The reference value of energy density is 1 × 10−12 J/m2.

Generate image description

Figure 19. The energy flow map of the energy distribution predicted by the method proposed in this article. Here, the frequency f = 487 Hz, and the damping η = 0.05. The reference value of energy flow is 1 × 10−12 W/m2.

Generate image description

Figure 20. The energy distribution along the diagonal of the plate. Here, the frequency f = 487 Hz and the damping η = 0.20. Exact (blue line) represents the solution given by the literature [42], EFEA (black line) represents the predicted value obtained from the energy finite element analysis, and Numerical (red line) represents the predicted value of the method proposed in this article. The reference value of energy density is 1 × 10−12 J/m2.

Generate image description

Figure 21. The energy distribution predicted by the method proposed in this article. Here, the frequency f = 487 Hz, and the damping η = 0.20. The reference value of energy density is 1 × 10−12 J/m2.

Generate image description

Figure 22. The energy flow map of the energy distribution predicted by the method proposed in this article. Here, the frequency f = 487 Hz, and the damping η = 0.20. The reference value of energy flow is 1 × 10−12 W/m2.

Table 1

Frequency and loss factor for each case study.

Case	Frequency (Hz)	Loss Factor
1	239	0.05
2	239	0.20
3	487	0.05
4	487	0.20

Appendix A. A Simple Example for Boundary Mapping

Consider the simplest square structure, as simple as only 4 elements and 12 half-boundaries. Make the order of polynomial and damping be 0, then the boundary mapping matrix T of will be a 12 × 12 sparse matrix. Table A1 shows the values of T.

Generate image description

Figure A1. Plate model (4 elements); red numbers represent half edges.

Table A1

Boundary mapping matrix T.

		Mapping Edge
		1	2	3	4	5	6	7	8	9	10	11	12
Receivingedge	1	0	0.5947	0.5947	0	0	0	0	0	0	0	0	0
	2	0	0	0	0.5947	0.2929	0	0	0	0	0	0	0
	3	0	0	0	0	0	0	0	0	0	0.5947	0	0.2929
	4	0	0	0	0	0.5947	0.5947	0	0	0	0	0	0
	5	0	0	0	0	0	0	0.5947	0.2929	0	0	0	0
	6	0.5947	0	0.2929	0	0	0	0	0	0	0	0	0
	7	0	0	0	0	0	0	0	0.5947	0.5947	0	0	0
	8	0	0	0	0	0	0	0	0	0	0.5947	0.2929	0
	9	0	0	0	0.5947	0	0.2929	0	0	0	0	0	0
	10	0	0	0	0	0	0	0	0	0	0	0.5947	0.5947
	11	0.5947	0.2929	0	0	0	0	0	0	0	0	0	0
	12	0	0	0	0	0	0	0.5947	0	0.2929	0	0	0

In fact, for any matrix, it represents a transformation law, and its leading eigenvalue represents the scaling ratio after this transformation. In this method, the boundary mapping matrix acts as a law to transform input energy to steady-state energy; in the case of 0 damping, the leading eigenvalue of T should be one (see Figure A2).

Generate image description

Figure A2. The eigenvalues of T. The circle represents 1, the eigenvalue that falls on the circle is 1, and the eigenvalue that falls within the circle is less than 1.

It can be seen that the leading eigenvalue falls on the unit circle. When the damping is not 0, the leading eigenvalue will fall in the unit circle. It should be pointed out that only when the damping is not 0 will the Neumann series meet the convergence condition, and at only this time will the steady-state solution exist. The leading eigenvalue can be understood simply as the ratio of the total energy of the system before and after mapping. When the system damping is 0, the leading eigenvalue is one, which means that there is no energy loss in the mapping process. When the system damping is not 0, the leading eigenvalue will be less than one, which means that energy is lost in the mapping process. In this method, when the system damping is 0, the Neumann series will not converge, which means that this method cannot simulate the vibration of the structure in the 0 damping state.

References

1. Yin, X.W.; Gui, H.F.; Gu, X.J.; Huang, F.; Shen, R.Y. Relevancy among power flow theory, statistical energy analysis and energy finite element method. J. Ship Mech.; 2007; 11, pp. 637-646.

2. Lyon, R.H.; Maidanik, G. Power flow between linearly coupled oscillators. J. Acoust. Soc. Am.; 1962; 34, pp. 623-639. [DOI: https://dx.doi.org/10.1121/1.1918177]

3. Lyon, R.H. Statistical analysis of power injection and response in structures and rooms. J. Acoust. Soc. Am.; 1969; 45, pp. 545-565. [DOI: https://dx.doi.org/10.1121/1.1911422]

4. Fahy, F.; De-Yuan, Y. Power flow between non-conservatively coupled oscillators. J. Sound Vib.; 1987; 114, pp. 1-11. [DOI: https://dx.doi.org/10.1016/S0022-460X(87)80227-4]

5. Beshara, M.; Keane, A. Statistical energy analysis of multiple, non-conservatively coupled systems. J. Sound Vib.; 1996; 198, pp. 95-122. [DOI: https://dx.doi.org/10.1006/jsvi.1996.0559]

6. Nefske, D.J.; Sung, S.H. Power flow finite element analysis of dynamic systems: Basic theory and application to beams. J. Vib. Acoust.; 1989; 111, pp. 94-100. [DOI: https://dx.doi.org/10.1115/1.3269830]

7. Glassner, A.S. An Introduction to Ray Tracing; Elsevier: Amsterdam, The Netherlands, 1989.

8. Kuttruff, H. Room Acoustics; CRC Press: Boca Raton, FL, USA, 2016.

9. McKown, J.; Hamilton, R. Ray tracing as a design tool for radio networks. IEEE Netw.; 1991; 5, pp. 27-30. [DOI: https://dx.doi.org/10.1109/65.103807]

10. Cerveny, V. Seismic Ray Theory; Cambridge University Press: Cambridge, UK, 2005.

11. Kulkarni, S.; Leppington, F.; Broadbent, E. Vibrations in several interconnected regions: A comparison of SEA, ray theory and numerical results. Wave Motion; 2001; 33, pp. 79-96. [DOI: https://dx.doi.org/10.1016/S0165-2125(00)00065-2]

12. Zhong, R.; Wang, Q.; Huang, Z.; Chen, L.; Shao, W.; Shuai, C. In-plane dynamic analysis of complex-shaped laminated cracked plates with irregular holes. AIAA J.; 2023; 61, pp. 3172-3189. [DOI: https://dx.doi.org/10.2514/1.J062774]

13. He, D.; Wang, Q.; Zhong, R.; Qin, B. Vibration analysis of functionally graded material (FGM) double layered floating raft structure by the spectro-geometric method. Structures; Elsevier: Amsterdam, The Netherlands, 2023; Volume 48, pp. 533-550.

14. He, D.; Wang, Q.; Zhong, R.; Qin, B. A unified analysis model of FGM double-layered submarine type coupled structure with spectral geometry method. Ocean Eng.; 2023; 267, 113213. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2022.113213]

15. Hu, S.; Zhong, R.; Wang, Q.; Qin, B.; Shao, W. A strong-form Chebyshev-RPIM meshless solution for free vibration of conical shell panels with variable thickness and fiber curvature. Compos. Struct.; 2022; 296, 115884. [DOI: https://dx.doi.org/10.1016/j.compstruct.2022.115884]

16. Zhong, S.; Jin, G.; Ye, T. Isogeometric vibration and material optimization of rotating in-plane functionally graded thin-shell blades with variable thickness. Thin-Walled Struct.; 2023; 185, 110593. [DOI: https://dx.doi.org/10.1016/j.tws.2023.110593]

17. Tian, L.; Jin, G.; He, T.; Ye, T.; Liu, Z.; Khadimallah, M.A.; Li, Z. A hybrid analytic–numerical formulation for the vibration analysis of a cylindrical shell coupled with an internal flexural floor structure. Thin-Walled Struct.; 2023; 183, 110382. [DOI: https://dx.doi.org/10.1016/j.tws.2022.110382]

18. Xue, Y.; Jin, G.; Zhang, C.; Han, X.; Chen, J. Free vibration analysis of functionally graded porous cylindrical panels and shells with porosity distributions along the thickness and length directions. Thin-Walled Struct.; 2023; 184, 110448. [DOI: https://dx.doi.org/10.1016/j.tws.2022.110448]

19. Chen, Y.; Ye, T.; Jin, G.; Lee, H.P.; Ma, X. A unified quasi-three-dimensional solution for vibration analysis of rotating pre-twisted laminated composite shell panels. Compos. Struct.; 2022; 282, 115072. [DOI: https://dx.doi.org/10.1016/j.compstruct.2021.115072]

20. Zhang, S.; Zhu, X.; Li, T.; Yin, C.; Li, Q.; Chen, R. A unified solution for free vibration analysis of beam-plate-shell combined structures with general boundary conditions. Int. J. Struct. Stab. Dyn.; 2022; 22, 2250080. [DOI: https://dx.doi.org/10.1142/S0219455422500808]

21. Zhang, S.; Li, T.; Zhu, X.; Yin, C.; Li, Q. Far field acoustic radiation and vibration analysis of combined shells submerged at finite depth from free surface. Ocean Eng.; 2022; 252, 111198. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2022.111198]

22. Nie, R.; Li, T.; Zhu, X.; Zhang, C. Analysis of free vibration characteristics of cylindrical shells with finite submerged depth based on energy variational principle. Symmetry; 2021; 13, 2162. [DOI: https://dx.doi.org/10.3390/sym13112162]

23. Nie, R.; Li, T.; Zhu, X.; Zhou, H. A general Fourier formulation for in-plane and out-of-plane vibration analysis of curved beams. Shock Vib.; 2021; 2021, 5511884. [DOI: https://dx.doi.org/10.1155/2021/5511884]

24. Xie, F.; Qu, Y.; Zhang, W.; Peng, Z.; Meng, G. Nonlinear aerothermoelastic analysis of composite laminated panels using a general higher-order shear de-formation zig-zag theory. Int. J. Mech. Sci.; 2019; 150, pp. 226-237. [DOI: https://dx.doi.org/10.1016/j.ijmecsci.2018.10.029]

25. Qu, Y.; Zhang, W.; Peng, Z.; Meng, G. Time-domain structural-acoustic analysis of composite plates subjected to moving dynamic loads. Compos. Struct.; 2019; 208, pp. 574-584. [DOI: https://dx.doi.org/10.1016/j.compstruct.2018.09.103]

26. Tornabene, F.; Fantuzzi, N.; Bacciocchi, M. Refined shear deformation theories for laminated composite arches and beams with variable thickness: Natural frequency analysis. Eng. Anal. Bound. Elements; 2019; 100, pp. 24-47. [DOI: https://dx.doi.org/10.1016/j.enganabound.2017.07.029]

27. Tornabene, F.; Viola, E. Free vibration analysis of functionally graded panels and shells of revolution. Meccanica; 2009; 44, pp. 255-281. [DOI: https://dx.doi.org/10.1007/s11012-008-9167-x]

28. Fantuzzi, N.; Tornabene, F.; Bacciocchi, M.; Dimitri, R. Free vibration analysis of arbitrarily shaped Functionally Graded Carbon Nanotube-reinforced plates. Compos. Part B Eng.; 2017; 115, pp. 384-408. [DOI: https://dx.doi.org/10.1016/j.compositesb.2016.09.021]

29. Żur, K.K.; Arefi, M.; Kim, J.; Reddy, J.N. Free vibration and buckling analyses of magneto-electro-elastic FGM nanoplates based on nonlocal modified higher-order sinusoidal shear deformation theory. Compos. Part B Eng.; 2019; 182, 107601. [DOI: https://dx.doi.org/10.1016/j.compositesb.2019.107601]

30. Mohammad-Rezaei Bidgoli, E.M.-R.; Arefi, M. Free vibration analysis of micro plate reinforced with functionally graded graphene nanoplatelets based on modified strain-gradient formulation. J. Sandw. Struct. Mater.; 2019; 23, pp. 436-472. [DOI: https://dx.doi.org/10.1177/1099636219839302]

31. Arefi, M.; Kiani, M.; Rabczuk, T. Application of nonlocal strain gradient theory to size dependent bending analysis of a sandwich porous nanoplate integrated with piezomagnetic face-sheets. Compos. Part B Eng.; 2019; 168, pp. 320-333. [DOI: https://dx.doi.org/10.1016/j.compositesb.2019.02.057]

32. Wu, Q.; Wang, S.; Yao, M.; Niu, Y.; Wang, C. Nonlinear dynamics of three-layer microplates: Simultaneous presence of the micro-scale and imperfect effects. Eur. Phys. J. Plus; 2024; 139, 446. [DOI: https://dx.doi.org/10.1140/epjp/s13360-024-05255-3]

33. Cen, Q.; Xing, Z.; Wang, Q.; Li, L.; Wang, Z.; Wu, Z.; Liu, L. Molding simulation of airfoil foam sandwich structure and interference optimization of foam-core. Chin. J. Aeronaut.; 2024; 37, pp. 325-338. [DOI: https://dx.doi.org/10.1016/j.cja.2024.08.025]

34. Wang, C.; Song, Z.; Fan, H. Novel evidence theory-based reliability analysis of functionally graded plate considering thermal stress behavior. Aerosp. Sci. Technol.; 2024; 146, 108936. [DOI: https://dx.doi.org/10.1016/j.ast.2024.108936]

35. Fu, T.; Hu, X.; Yang, C. Impact response analysis of stiffened sandwich functionally graded porous materials doubly-curved shell with re-entrant honeycomb auxetic core. Appl. Math. Model.; 2023; 124, pp. 553-575. [DOI: https://dx.doi.org/10.1016/j.apm.2023.08.024]

36. Yu, Y.; Fu, T.; Wang, S.; Yang, C. Dynamic response of novel sandwich structures with 3D sinusoid-parallel-hybrid honeycomb auxetic cores: The cores based on negative Poisson’s ratio of elastic jump. Eur. J. Mech.-A/Solids; 2025; 109, 105449. [DOI: https://dx.doi.org/10.1016/j.euromechsol.2024.105449]

37. Li, H.; Xu, J.; Pang, F.; Gao, C.; Zheng, J. A unified Jacobi-Ritz-spectral BEM for vibro-acoustic behavior of spherical shell. Comput. Math. Appl.; 2024; 176, pp. 415-431. [DOI: https://dx.doi.org/10.1016/j.camwa.2024.10.031]

38. Li, H.; Xu, J.; Gong, Q.; Teng, Y.; Pang, F.; Zhang, L. A study on the dynamic characteristics of the stiffened coupled plate with the effect of the dynamic vibration absorbers. Comput. Math. Appl.; 2024; 168, pp. 120-132. [DOI: https://dx.doi.org/10.1016/j.camwa.2024.04.026]

39. Pang, F.; Tang, Y.; Li, C.; Zhou, F.; Gao, C.; Li, H. Reconstructed source method for underwater noise prediction of a stiffened cylindrical shell. Ocean Eng.; 2024; 310, 118828. [DOI: https://dx.doi.org/10.1016/j.oceaneng.2024.118828]

40. Cvitanovic, P.; Artuso, R.; Mainieri, R.; Tanner, G.; Vattay, G.; Whelan, N.; Wirzba, A. Chaos: Classical and quantum. ChaosBook. Org; Niels Bohr Institute: Copenhagen, Denmark, 2005; Volume 69, 25.

41. Meyer, C.D. Matrix Analysis and Applied Linear Algebra; Society for Industrial and Applied Mathematics: Alexandria, VA, USA, 2000; Volume 71.

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