Dongze Cui¹, Rafael Da Silva Raqueti², Najib Fazail³, Da Chen^1
1Tianjin Engineering Research Center of Civil Aviation Energy Environment and Green Development, Civil Aviation University of China, Tianjin, China.
²Groupe d’Acoustigue de l’ Université de Sherbrooke (GAUS), Université de Sherbrooke, de Université, Sherbrooke, Québec, Canada.
³Vibratec Asia Pacifc, Jalan Bangsar Utama 1, Kuala Lumpur, Malaysia.
DOI: 10.4236/jamp.2025.138162   PDF    HTML   XML   95 Downloads   417 Views   Citations

Abstract

A typical inverse problem, identifying the wavespace for metastructures from measurements of displacement fields, is crucial to understanding the dynamic properties of metastructures and the optimal design of vibroacoustic performance. The validity examination of inverse methods for heterogeneous metastructures is still an open question, especially for those showing wave coupling effects. The present study examines the robustness and accuracy of different inverse methods for extracting the complex wavespace. We compare the nonlinear Inhomogeneous Wave Correlation (IWC) method and its extension, the Green’s Function Correlation (GFC) framework, with the linear Algebraic Wavenumber Identification (AWI) framework. Direct wave-based methodologies are employed to assess both the precision of these techniques and the influence of wave coupling on the inverse methods’ robustness. The study starts by evaluating the precision of all methods through Highly Contrasted Dissipative Structures (HCDS) with rheological viscoelastic material. Then, the limitation of GFC due to the assumption of the elastic formulation is demonstrated by the sandwich plate with thick soft cores.

Keywords

Heterogeneous Metastructures, Wave Propagation, Inverse Wavespace Identification Frameworks, Green’s Function Correlation (GFC), Algebraic Wavenumber Identification (AWI), Damping Loss Factor (DLF)

Share and Cite:

Keywords:

1. Introduction

Heterogeneous metastructures are engineering structures designed and manufactured artificially with vibroacoustic performance that significantly surpasses traditional structures. Metastructures are widely used in the transport engineering sector, such as the mainframes of automobiles, aircraft structures, high-speed railways, and naval architectures, commonly including sandwich composite structures, layered composite materials, and stiffened plates. By precisely controlling the geometry, size, material distribution, and arrangement of periodic units, these structures achieve the suppression and manipulation of wave propagation (such as elastic waves and acoustic waves) within specific frequency ranges [1]-[5].

A complete closed-loop vibroacoustic design framework of metastructure involves two complementary approaches: the direct design method and the reverse validation approach. Each path serves a distinct purpose in ensuring the performance and functionality of the heterogeneous metastructures.

The direct design method starts from material properties and geometric layout of the metastructure to predict its vibroacoustic indicators through analytical or numerical approaches. Typical indicators include dispersion relations, and Damping Loss Factor (DLF). The direct design method then optimizes the vibroacoustic performance of the metastructure by adjusting structural and material properties to meet the vibration and noise reduction requirements within a specific frequency range.

Common methods include the Finite Element Method (FEM), a versatile numerical technique used to model and analyze the physical behavior of metastructures under various loading conditions [6]-[14]. FEM discretizes the fluid, solid, and acoustic fields and solves differential equations based on conservation laws of momentum, mass, or energy under given boundary conditions [15], and approximates field variables using local, typically polynomial, predefined shape functions, providing accurate estimations of acoustic indicators [16]. FEM can accurately simulate the interlayer multi-scale dynamic behavior of laminated structures and accurately calculate energy transmission in metastructures. However, FEM involves extensive computations, and as frequency increases, the accuracy of its approximate shape functions and model complexity reduces computational efficiency [15] [17]. FEM also encounters challenges due to approximate shape functions and increasing model complexity with frequency [18]-[20]. To enhance computational efficiency, formulations proposed in [21]-[23] are under consideration. An efficient alternative is the Wave-based FEM (WFEM), which uses only the Representative Unit Cell (RUC) of the metastructures and employs the Bloch-Floquet theorem to simulate the boundary conditions to efficiently compute the vibroacoustic indicators [1]-[5].

Contrary to the direct design methods, wavespace identification is a classical inverse problem. It identifies the dispersion curves from the structural forced response of the metastructure via mathematical relations between spatial response spectrum into the wavenumber-frequency domain. This process further identifies the bandgaps and equivalent structural parameters, such as elastic modulus and DLF. The wavespace identification technique offers an inverse validation approach for vibration and noise reduction design in metastructures, enabling closed-loop design within the Bloch wave theory framework. This approach enhances the efficiency and reliability of metastructure design. Furthermore, by identifying the wave propagation characteristics within metastructures, the inverse approach complements direct design methods, elucidating the dynamic behavior and wave transmission mechanisms of complex structures.

For inverse approaches, structural responses under random or harmonic excitations obtained from experimental or numerical methods can be used to estimate the complex wavenumber space, which contains information on the dispersion relations, energy propagation, and modal properties of metastructures, to propose novel DLF estimation methods [24] [25]. Waveapace identification frameworks are consequently developed to predict the dispersion relations and DLF of metastructures. Existing methods predominantly focus on the real part of the wavenumber, resulting in an inadequate prediction of the DLF, which is linked to its imaginary part. Two generic approaches corresponding to the linear and nonlinear approaches are the Algebraic Wavenumber Identification (AWI) method based on the Laplace transforms [26], and the Inhomogeneous Wave Correlation (IWC) method based on the Fourier transform [24] [25] [27].

The classical methods for DLF estimation include the modal technique (also known as the half-power bandwidth method) [28], the Decay Rate Method (DRM) [29], and the Power Input Method (PIM) [15], which is based on the principle of energy conservation by equating the input power and dissipated energy [30]. PIM has advantages over other methods, such as being independent of mode shapes or natural frequencies and accommodating multiple modes and non-linearities [24].

The IWC method inverts the wave propagation characteristics of a two dimensional structures from FRF measurements by maximizing the correlation between an inhomogeneous wave and the displacement field. However, the plane wave assumption in IWC is invalid near the excitation point due to evanescent waves. Therefore, the measurement window must be far from the loaded zone, which is impractical in real-world scenarios. To overcome this limitation, Tufano et al. [27] proposed the Green’s Function Correlation (GFC) method, which uses a Green’s function-based model to generate an inhomogeneous wave that accounts for the structure properties. The GFC method was applied to an isotropic laminated plate and an isotropic plate with tuned mass dampers. However, both IWC and GFC methods require solving a nonlinear wavenumber search problem, which is computationally expensive.

The AWI technique is a linear method that can extract the complex wavenumbers of wave propagation in periodic structures from FRF obtained experimentally or numerically. The AWI technique is based on the algebraic approach of parameter identification: The algebraic derivatives method is applied to the spatial displacement field of the structure at each frequency. The Laplace transform is used to convert the differential equation that governs the wave propagation into an algebraic equation. The inverse Laplace transform converts the algebraic equation back into the spatial domain, resulting in a new linear regression equation with multiple integrals. The complex wavenumbers are estimated by solving the multiple integrals, thereby improving the computational efficiency of AWI. The AWI technique is computationally efficient and robust to noise and measurement errors, and can be applied to planar and cylindrical periodic structures [26].

In this paper, firstly, nonlinear and linear wavespace identification frameworks are presented to identify the complex wavenumbers from the displacement field obtained by Finite Element Method (FEM). Secondly, a DLF estimation method using complex wavenumbers is developed, which also derives the average DLF of non-isotropic metastructures using the modal density. Finally, numerical examples of various metastructures, especially the sandwich structures with inhomogeneous cores, are validated by the Wave Finite Element (WFE) scheme to assess the accuracy and efficiency of the proposed method.

The paper is organized as follows: Section 2 introduces the nonlinear and linear inverse problem algorithms for retrieving the complex wavenumbers from the FEM displacement field. Section 3 examines the presented methods with various planar structures, including a sandwich structure with a thick soft core, and a highly contrasted and dissipative metastructure with viscoelastic cores. Section 4 discusses and concludes the results.

2. Overview of Nonlinear and Linear Wave Identification Methods

In practical applications, acquisition grids are generally distributed uniformly over the two-dimensional surface of the metastructure. Displacement measurements at these 2D grid points are taken at specific angles to facilitate wave inversion in the intended propagation directions, as illustrated in Figure 1. Using inverse problem algorithms, the complex wavenumber $\kappa \left(\omega ,\theta \right)$ is estimated, allowing the extraction of the Damping Loss Factor (DLF).

2.1. Nonlinear Wavespace Identification Framework

2.1.1. Inhomogeneous Wave Correlation

The classical two-dimensional IWC aims to choose an inhomogeneous plane wave in the polar coordinate system. The equation describing the properties of the wavenumber, $\kappa$ , propagating in direction $\theta$ is defined as:

${\hat{\sigma }}_{\kappa ,\gamma ,\theta }\left(x_i,y_j\right)=e^{-i\kappa \left(\theta \right)\left(1+i\gamma \left(\theta \right)\right)\left(x_{i}cos\theta +y_{i}sin\theta \right)}$ (1)

where $\gamma$ denotes the attenuation factor, the complex wavenumber can also be expressed as $\hat{\kappa }\left(\theta \right)=\kappa \left(\omega ,\theta \right)\left(1+i\gamma \left(\omega ,\theta \right)\right)$ , $\left(x_i,y_j\right)$ is the coordinate of a random acquisition point.

The IWC is based on searching for the maximum of the correlation function between the measured displacement field $\hat{w}\left(x,y\right)$ and the function parameterized by the complex wave number. The correlation function on the spatial domain Ω

writes [31]:

$\text{IWC}\left(\kappa ,\gamma ,\theta \right)=\frac{\left|{\displaystyle {\iint }_{\Omega }\hat{w}\left(x_i,y_i\right)\cdot {\hat{\sigma }}_{\kappa ,\gamma ,\theta }^{*}\left(x_i,y_i\right)dxdy}\right|}{\sqrt{{\displaystyle {\iint }_{\Omega }{\left|\hat{w}\left(x_i,y_i\right)\right|}^{2}dxdy}\cdot {\displaystyle {\iint }_{\Omega }{\left|{\hat{\sigma }}_{\kappa ,\gamma ,\theta }\left(x_i,y_i\right)\right|}^{2}dxdy}}}$ (2)

where $\ast$ denotes the complex conjugate.

The double integrations must be approximated numerically to proceed with the measured displacement field on a discrete subset of Ω. Rewriting Equation (2) in the discrete domain, the integration over the entire surface Ω is replaced by a finite weighted sum:

${\displaystyle {\iint }_{\Omega }dxdy}=\underset{j=1}{\overset{N}{{\displaystyle \sum }}}{\rho }_j{\Omega }_j$ (3)

where ${\rho }_j$ is the coherence of the measured signal at each point ( ${\rho }_j=1$ if the coherence is not available), ${\Omega }_j$ is an estimation of the surface around the point $j$ and $N$ is the total number of acquisition points.

When measurement grids are known, it is preferable to incorporate the grid information into the numerical approximation of the scalar product and the norm integrals [32] [33]. The coherence of the measurement associated with the $\left(i,j\right)$ -th grid surface $S_{i,j}$ is denoted as ${\rho }_{i,j}$ , Equation (2) then becomes:

$\text{IWC}\left(\kappa ,\gamma ,\theta \right)=\frac{\left|{{\displaystyle \sum }}_{i=1}^N{{\displaystyle \sum }}_{j=1}^M\hat{w}\left(x_i,y_j\right){\hat{\sigma }}_{\kappa ,\gamma ,\theta }^{*}\left(x_i,y_j\right){\rho }_i{S}_i\right|}{\sqrt{{{\displaystyle \sum }}_{i=1}^N{{\displaystyle \sum }}_{j=1}^M{\left|\hat{w}\left(x_i,y_j\right)\right|}^2{\rho }_i{S}_i{{\displaystyle \sum }}_{i=1}^N{{\displaystyle \sum }}_{j=1}^M{\left|{\hat{\sigma }}_{\kappa ,\gamma ,\theta }\left(x_i,y_j\right)\right|}^2{\rho }_i{S}_i}}$ (4)

where ${\rho }_i$ is regarded as the surface integration weight at $i$ -th grid, $S_i$ is an estimation of the surface around the point $i$ and $N$ is the total number of acquisition points [32] [33].

2.1.2. Correlation Model with Green’s Function

This section defines an improved correlation model: the inhomogeneous wave model is replaced by Green’s function-based model to simulate the dynamic behavior of an infinite Kirchhoff-Love plate. The Green’s function for the measured displacement field on a thin plate of infinite dimensions is given by [34] [35]:

The equation describing the properties of the wavenumber, $\kappa$ , propagating in direction $\theta$ is defined as:

${\hat{\sigma }}_{\kappa ,\gamma }=G_{\infty }\left(\hat{\kappa },r\right)=\frac{i}{8{\hat{\kappa }}^{2}D}\left[H_0^{\left(1\right)}\left(\hat{\kappa }{r}_i\right)-H_0^{\left(1\right)}\left(i\hat{\kappa }{r}_i\right)\right]$ (5)

where $G_{\infty }$ denotes the Green’s function associated with an infinite plate, the complex wavenumber $\hat{\kappa }$ is defined as $\hat{\kappa }={\kappa }_R+i{\kappa }_I=\kappa \left(1+\gamma \right)$ , the flexural stiffness is defined as $D=\frac{E{h}^3}{12\left(1-{\nu }^2\right)}$ , with $E$ representing Young’s modulus, $h$ designating the thickness, $\nu$ representing Poisson’s coefficient, radius $r_i=\sqrt{{(x_i-x_e)}^2+{(y_i-y_e)}^2}$ is defined as the spatial distance between the excitation point $\left(x_e,y_e\right)$ and the observation point $\left(x_i,y_i\right)$ .

In this case, the correlation function to be maximized writes:

$\text{GFC}\left(\hat{\kappa },r,\theta \right)=\frac{\left|{\displaystyle {\iint }_{\Omega }\hat{w}\left(r,\theta \right)\cdot {\hat{\sigma }}_{\kappa ,\gamma }^{*}\left(\hat{\kappa },r,\theta \right)dxdy}\right|}{\sqrt{{\displaystyle {\iint }_{\Omega }{\left|\hat{w}\left(r,\theta \right)\right|}^{2}dxdy}\cdot {\displaystyle {\iint }_{\Omega }{\left|{\hat{\sigma }}_{\kappa ,\gamma }\left(\hat{\kappa },r,\theta \right)\right|}^{2}dxdy}}}$ (6)

To facilitate the following analysis, it is preferable to eliminate the contribution of the flexural stiffness, as expressed in Equation (5), by introducing the dispersion relation about the flexural wavenumber using the Kirchhoff-Love thin plate theory:

${\hat{\kappa }}^2=\sqrt{\frac{m_s}{D}}\omega$ (7)

where $m_s$ is the mass per unit area and $\omega$ is the angular frequency. The GFC model becomes [25]:

$G\left(\hat{\kappa },r,\theta \right)=\frac{i{\hat{\kappa }}^2\left(\theta \right)}{8m_s{\omega }^2}\left[H_0^{\left(1\right)}\left(\hat{\kappa }\left(\theta \right)r_i\right)-H_0^{\left(1\right)}\left(i\hat{\kappa }\left(\theta \right)r_i\right)\right]$ (8)

The function offers a means to evaluate the equivalent elastic properties of intricate structures under various propagation angles.

In polar coordinates, the GFC model is expressed as follows:

$\text{GFC}\left(\hat{\kappa },r,\theta \right)=\frac{\left|{{\displaystyle \sum }}_{i=1}^N\hat{w}\left(r_i,\theta \right)G^{*}\left(\hat{\kappa },r_i,\theta \right){\rho }_i{S}_i\right|}{\sqrt{{{\displaystyle \sum }}_{i=1}^N{\left|\hat{w}\left(r_i,\theta \right)\right|}^2{\rho }_i{S}_i{{\displaystyle \sum }}_{i=1}^N{\left|G\left(\hat{\kappa },r_i,\theta \right)\right|}^2{\rho }_i{S}_i}}$ (9)

Note that the same simplification from double integration to summation is inherently conducted in the polar coordinate system. The determination of the complex wavenumber is achieved by maximizing the function $\text{GFC}\left(\hat{\kappa },r,\theta \right)$ for each angle and frequency.

2.2. Linear Wavespace Identification Framework

2.2.1. Algebraic Wavenumber Identification

The harmonic displacement at any measurement point $\left(x_n,y_n\right)$ on a plate can be effectively modeled through the superposition of $n_{\omega }$ plane waves. Given the wave propagation angle $\theta$ , the displacement at any measurement point in the wave propagation direction can be expressed as follows:

$\hat{w}\left(\theta ,{\hat{r}}_n\right)=\underset{m=1}{\overset{n_{\omega }}{{\displaystyle \sum }}}{A}_m{e}^{p_{\theta ,m}{\hat{r}}_n}=\underset{m=1}{\overset{n_{\omega }}{{\displaystyle \sum }}}{A}_m{e}^{i{\kappa }_{\theta ,m}{\hat{r}}_n}$ (10)

In the wavenumber domain, the Laplace transform of the displacement field can be expressed as follows:

$W\left(\theta ,s\right)=\underset{m=1}{\overset{n_{\omega }}{{\displaystyle \sum }}}\frac{A_m}{s-p_{\theta ,m}}$ (11)

The measured displacement in polar coordinates can be viewed as a solution to an Ordinary Differential Equation (ODE). In the wavenumber domain, the characteristic polynomial of this ODE can be expressed as follows:

$\Psi \left(s\right)=\underset{m=1}{\overset{n_{\omega }}{{\displaystyle \prod }}}\left(s-p_{\theta ,m}\right)=\underset{i=0}{\overset{n_{\omega }}{{\displaystyle \sum }}}\gamma \left(n_{\omega }-i\right)s^i$ (12)

where $\gamma {\left(i\right)}_{i\in \left[0,n_{\omega }\right]}$ represent the unknown coefficients of the characteristic polynomial. The wavenumber can be determined by solving the polynomial provided all the coefficients. Consequently, a novel function in the wavenumber domain is formulated by the multiplication of Equations (11) and (12):

$W\left(\theta ,s\right)\Psi \left(s\right)=\underset{m=1}{\overset{n_{\omega }}{{\displaystyle \sum }}}\frac{A_m}{s-p_{\theta ,m}}\underset{m=1}{\overset{n_{\omega }}{{\displaystyle \prod }}}\left(s-p_{\theta ,m}\right)=\underset{m=1}{\overset{n_{\omega }}{{\displaystyle \sum }}}{A}_m\underset{i=1,i\ne m}{\overset{n_{\omega }}{{\displaystyle \prod }}}s-p_{\theta ,i}$ (13)

The $n_{\omega }$ -th differential equations of $\left(n_{\omega }-1\right)$ -th order polynomial $U\left(\theta ,s\right)\Psi \left(s\right)$ with respect to $s$ writes:

$\frac{d^{n_{\omega }}\left[W\left(\theta ,s\right)\Psi \left(s\right)\right]}{d{s}^{n_{\omega }}}=\frac{d^{n_{\omega }}\left[W\left(\theta ,s\right){{\displaystyle \sum }}_{i=0}^{n_{\omega }}\gamma \left(n_{\omega }-i\right)s^i\right]}{d{s}^{n_{\omega }}}=0$ (14)

To compute this equation, the Leibniz formula is employed:

$\frac{d^{n_{\omega }}\left[W\left(\theta ,s\right)\Psi \left(s\right)\right]}{d{s}^{n_{\omega }}}=\underset{j=0}{\overset{n_{\omega }}{{\displaystyle \sum }}}\left(\begin{array}{c}{n}_{\omega }\\ j\end{array}\right)\frac{d^{n_{\omega }-j}W\left(\theta ,s\right)}{d{s}^{n_{\omega }-j}}\frac{d^j\Psi \left(s\right)}{d{s}^j}$ (15)

and the high-order algebraic derivatives are defined as follows:

$\frac{d^{n_{\omega }}\left(s^j\right)}{d{s}^{n_{\omega }}}=\frac{j!}{\left(j-n_{\omega }\right)!}{s}^{j-n_{\omega }}$ (16)

where $!$ symbolizes the factorial operation of an integer.

The ensuing equation can be derived as follows:

$\frac{d^{n_{\omega }}\left(s^{i}W\left(\theta ,s\right)\right)}{d{s}^{n_{\omega }}}=s^{i-n_{\omega }}\underset{j=0}{\overset{n_{\omega }}{{\displaystyle \sum }}}\frac{d^{j}W\left(\theta ,s\right)}{d{s}^j}{s}^i\left(\begin{array}{c}i\\ n_{\omega }-j\end{array}\right)\left(n_{\omega }-j\right)!\left(\begin{array}{c}{n}_{\omega }\\ j\end{array}\right)$ (17)

Integrate the provided equation with the $n_{\omega }$ -th differential equation as follows:

$\underset{i=0}{\overset{n_{\omega }}{{\displaystyle \sum }}}\underset{j=i}{\overset{n_{\omega }}{{\displaystyle \sum }}}\left(\begin{array}{c}{n}_{\omega }\\ j\end{array}\right)\left(\begin{array}{c}{n}_{\omega }-i\\ n_{\omega }-j\end{array}\right)\left(n_{\omega }-j\right)!s^{j-i}\frac{d^{j}W\left(\theta ,s\right)}{d{s}^j}\gamma \left(i\right)=0$ (18)

Subsequently, the $n_{\omega }$ -th differential equation is transformed back into the spatial domain using the Inverse Laplace Transform:

${\mathcal{L}}^{-1}\left(\frac{1}{s^I}\frac{d^{J}W\left(\theta ,s\right)}{d{s}^J}\right)=\frac{1}{\left(I-1\right)!}{\displaystyle {\int }_0^{r_n}{v}_{I-1,J}\left(\tau \right)W\left(\theta ,\tau \right)d\tau }$ (19)

with

$v_{I,J}\left(\tau \right)={\left(r_n-\tau \right)}^I{\left(-\tau \right)}^J$ (20)

In order to ensure compatibility with the Inverse Laplace Transform, the $n_{\omega }$ -th differential equation is divided by $s^{n_{\omega }+1}$ :

$\underset{i=0}{\overset{n_{\omega }}{{\displaystyle \sum }}}\underset{j=i}{\overset{n_{\omega }}{{\displaystyle \sum }}}\left(\begin{array}{c}{n}_{\omega }\\ j\end{array}\right)\left(\begin{array}{c}{n}_{\omega }-i\\ n_{\omega }-j\end{array}\right)\left(n_{\omega }-j\right)!\frac{1}{s^{n_{\omega }+1+i-j}}\frac{d^{j}U\left(\theta ,s\right)}{d{s}^j}\gamma \left(i\right)=0$ (21)

Upon the application of the Inverse Laplace Transform, the equation in the spatial domain is expressed as follows:

$\begin{array}{c}0=\underset{i=0}{\overset{n_{\omega }}{{\displaystyle \sum }}}\underset{j=i}{\overset{n_{\omega }}{{\displaystyle \sum }}}\left(\begin{array}{c}{n}_{\omega }\\ j\end{array}\right)\left(\begin{array}{c}{n}_{\omega }-i\\ n_{\omega }-j\end{array}\right)\left(n_{\omega }-j\right)!\frac{1}{\left(n_{\omega }+i-j\right)!}\\ \times {\displaystyle {\int }_0^{r_n}{\left(r_n-\tau \right)}^{n_{\omega }+i-j}{\left(-\tau \right)}^{j}U\left(\theta ,\tau \right)d\tau }\gamma \left(i\right)\end{array}$ (22)

It can be rewritten as follows:

$\underset{i=0}{\overset{n_{\omega }}{{\displaystyle \sum }}}\varphi \left(i,\theta ,r_n\right)\gamma \left(i\right)=0$ (23)

with

$\begin{array}{c}\varphi \left(i,\theta ,r_n\right)=\underset{i=0}{\overset{n_{\omega }}{{\displaystyle \sum }}}\underset{j=i}{\overset{n_{\omega }}{{\displaystyle \sum }}}\left(\begin{array}{c}{n}_{\omega }\\ j\end{array}\right)\left(\begin{array}{c}{n}_{\omega }-i\\ n_{\omega }-j\end{array}\right)\left(n_{\omega }-j\right)!\frac{1}{\left(n_{\omega }+i-j\right)!}\\ \times {\displaystyle {\int }_0^{r_n}{\left(r_n-\tau \right)}^{n_{\omega }+i-j}{\left(-\tau \right)}^{j}U\left(\theta ,\tau \right)d\tau }\end{array}$ (24)

where the numerical integration can be performed using the trapezoidal rule.

The third step involves the estimation of $\gamma \left(i\right)$ using the Least Squares method:

$HR=\left[\begin{array}{cccc}\varphi \left(n_{\omega },\theta ,r_1\right)& \varphi \left(n_{\omega }-1,\theta ,r_1\right)& \cdots & \varphi \left(0,\theta ,r_1\right)\\ \varphi \left(n_{\omega },\theta ,r_2\right)& \varphi \left(n_{\omega }-1,\theta ,r_2\right)& \cdots & \varphi \left(0,\theta ,r_2\right)\\ ⋮& ⋮& \ddots & ⋮\\ \varphi \left(n_{\omega },\theta ,r_n\right)& \varphi \left(n_{\omega }-1,\theta ,r_n\right)& \cdots & \varphi \left(0,\theta ,r_n\right)\end{array}\right]\left[\begin{array}{c}\gamma \left(n_{\omega }\right)\\ \gamma \left(n_{\omega }-1\right)\\ ⋮\\ \gamma \left(0\right)\end{array}\right]=0$ (25)

where $R$ represents the eigenvector associated with the smallest eigenvalue of the convolution of matrices $H^{*}H$ .

Upon acquiring the coefficient vector $R$ for the characteristic polynomial, the wavenumber at the propagation direction $\theta$ is determined as ${\kappa }_{\theta ,m}=-i{p}_{\theta ,m}$ .

2.2.2. Wavenumber Filter by Physical Constraints

Conventional wave filtering technique relies on the sign of the real part and the ratio of the imaginary to the real part, such technique easily loses robustness when more than four solutions are present. It is therefore strongly recommended to filter the wave output from the linear AWI approach by enforcing the continuity of group and phase velocities in the frequency domain:

$C_g=\frac{\partial \omega }{\partial \kappa },C_{\phi }=\frac{\omega }{\kappa }$ (26)

where group velocity $C_g$ signifies the speed of energy transmission of the wave in structures, phase velocity $C_{\phi }$ denotes the speed and direction at which the phase of a wave propagates through space [36].

The wavenumber can be filtered by a dynamic programming algorithm with physical constraints [37], with physical constraints defined by the continuity of group/phase velocities in the frequency domain, which are also computed by the wavenumbers.

The aim of the algorithm is to find the shortest path problem with physical constraints. The problem involves identifying a path from low to high frequency in a scatter plot, ensuring smooth transitions in group and phase velocities to maintain physical continuity. Dynamic programming is well-suited for this task, as it optimizes multi-stage decisions through state transitions and cost minimization, making it effective for mode tracking in noisy data.

The workflow of the global dynamic programming algorithm for modal tracking is summarized as follows:

1. Input and Initialization

The algorithm accepts a frequency vector ( $N\times 1$ ), a candidate wavenumber matrix ( $N\times M$ ), and a user-specified initial wavenumber. A dynamic programming cost matrix is initialized and a backtracking pointer matrix (parent) is set up. The candidate at the first frequency point closest to the initial wavenumber is assigned a zero cost, while all other candidates remain unreachable (cost = $\infty$ ).

2. Dynamic Programming Loop

For each consecutive frequency step (from $i=1$ to $N-1$ ), the algorithm computes the cost of transitioning from each candidate at frequency i to every candidate at frequency $i+1$ . The cost function is defined as:

${\text{cost}}_g=\alpha \left|v_i-v_{i+1}\left|+\left(1-\alpha \right)\right|{\kappa }_i-{\kappa }_{i+1}\right|$ (27)

where $v$ can be chosen as $C_g$ , $C_{\phi }$ or both, $\alpha$ is a weighting factor between 0 and 1. The cumulative cost is updated, and the corresponding parent pointer is recorded if a lower cost path is found.

3. Final Selection and Backtracking

At the final frequency point, the candidate with the minimal cumulative cost is selected. The optimal continuous modal branch is then retrieved by backtracking through the parent pointers from the final frequency point to the first.

This workflow above ensures that the filtered wave is physically continuous in the frequency range, mitigating local mismatches by considering global optimality.

3. Results of Homogeneous Metastructures

This section introduces the application of the proposed nonlinear and linear techniques to several numerical examples of varying complexity. Simple isotropic plate and orthotropic plate where analytical results are available. For heterogeneous meta-structures such as the sandwich laminate with a thick core and the orthotropic graphite-epoxy sandwich with a thin core, validation is conducted using the reference Wave Finite Element (WFE) scheme, and the numerical PIM based on flexural response data obtained from the full FEM analysis [15].

The computations were performed using MATLAB R2023a on an ASUS computer equipped with an Intel® Core^TM i7-10875H CPU @ 2.30 GHz and 16.0 GB of RAM, supplemented by two additional 4TB SSDs to enhance computing memory capacity. Furthermore, the computational efficiency of the WFE method is compared against FEM.

The material properties used for simulating the homogeneous metastructures are listed in Figure 1.

3.1. Sandwich Plate with a Thick Soft Core: The Impact of Symmetric Motion

A specific HCS instance is presented with a configuration of 2 mm skins and a 20 mm core [38]. The WFE scheme’s performance for soft thick cores is exemplified through this configuration for the examination of symmetric and asymmetric motions of HCS. The material properties are detailed in Table 1.

The angular frequency corresponding to the symmetric motion can be estimated analytically:

${\omega }_{sym}=\sqrt{\frac{E_s\ast \left(1-{\nu }_c\right)}{\left(1+{\nu }_c\right)\left(1-2{\nu }_c\right)}\frac{\frac{1}{m_{s,1}}+\frac{1}{m_{s,2}}}{h_c}}$ (28)

where subscripts $s$ and $c$ correspond to the skin and core, $m$ denotes the surface density, $h$ denotes the thickness.

As depicted in Figure 2, the wavenumber filter works robustly with the physical constraints. Besides, AWI works well for the inverse identification of flexural wavenumbers from the structural responses, while GFC outputs show a high discrepancy compared to the reference results.

The real and imaginary parts of the wavenumber linked to asymmetric and symmetric motions are depicted in Figure 3. Below the first out-of-phase symmetric mode at frequency $f_{sym}=281.4955$ Hz calculated using Equation (28), the wavenumber corresponding to the symmetric motion is evanescent with real and imaginary parts of similar magnitude, above this frequency, the wavenumber of the symmetric wavemode becomes propagative as its imaginary part shifts to zero.

The asymmetric and symmetric wavemodes are depicted in Figure 3 for a visible comprehension of the dilatational/compressional modeshape, the length of the FE model equals the wavelength of the target wave. The associated wavemode is asymmetric: at low frequency, the flexural motion is purely asymmetric, at high frequency, the wavemode is dominated by asymmetric flexural wave with transverse shear [39]-[41]. Therefore, consideration of the shear effect in the asymmetric motion is needed for this specific configuration family to compute the vibroacoustic indicators.

GFC doesn’t work for this structure due to its presumption of the elastic formulation inherently presented in Equation (7), which soon fails for metastructures showing dilatational symmetric motion in the out-of-plane direction [4].

3.2. Highly Contrasted and Dissipative Metastructures with a Rheological Core

The complex elastic modulus of SMP is governed by a relationship proposed and substantiated by Butaud [42], which writes as follows:

$E\left(\omega ,T_0\right)=E_0+\frac{\left(E_{\infty }-E_0\right)}{\left(1+\gamma {\left(i\omega \tau \right)}^{-k}+{\left(i\omega \tau \right)}^{-h_i}+{\left(i\omega \beta \tau \right)}^{-1}\right)}$ (29)

where $E_0=0.67$ , $E_{\infty }=2211$ , $k=0.16$ , $h_i=0.79$ , $\gamma =1.68$ , $\beta =3.8\times {10}^4$ , and ${\tau }_0=0.61$ .

Furthermore, the behavior of SMP adheres to the concept of time-temperature superposition, as observed in many polymers [43]. This phenomenon involves the characteristic time $t_T$ , which is connected to $t_{T_0}$ (the characteristic time at the reference temperature $T_0$ ) through a shifter $a_T$ governed by the following relation:

$log\left(a_T\right)=-C_1\frac{T-T_0}{C_2+T-T_0}$ (30)

where $C_1=10.87˚\text{C}$ , $C_2=32.57˚\text{C}$ , and the reference temperature $T_0=40˚\text{C}$ .

The weak elasticity and high damping ratio of SMP introduce complexity in numerical computations [42], the attenuation characteristics of plane waves are beneficial for testing the validity of WIM for their accuracy.

A sandwich structure with a thick tBA/PEGDMA core (more commonly referred to as SMP) core is simulated by the in-house FE package. The configuration consists of 0.5 mm aluminum skins, while the SMP core has a thickness of 2.2 mm, the material properties are listed in Table 1.

3.2.1. Shape Memory Polymer at 65˚C

Figure 4 displays the bending wavenumbers and DLF for the SMP65˚C sandwich plate. Notably, the bending wavenumbers computed using various methods exhibit excellent agreement, further substantiating the accuracy and reliability of the analyses.

The DLF values from the WIM are depicted alongside those from the AHM and PIM-FEM approaches [2]. This agreement of outcomes across different techniques supports the efficacy of the two WIMs in capturing the dispersion relation and DLF of the SMP65˚C sandwich panel.

The DLF displays excellent agreement with both the AHM scheme and the PIM-FEM outcomes [15]. It’s worth noting that discrepancies arising in the

low-frequency domain in the GFC results have been attributed to the inherent nature of the WIM method, which relies on the presence of multiple wavelengths in the out-of-plane displacement field to capture the wave propagation effectively. Moreover, the GFC approach inherently neglects the modal behavior of the panel in the free field Green’s function that is employed within the methodology [24].

3.2.2. Shape Memory Polymer at 80˚C

Figure 5 displays the bending wavenumbers and DLF for the SMP80˚C sandwich plate. Notably, the bending wavenumbers computed using various methods exhibit excellent agreement, further substantiating the accuracy and reliability of the analyses.

4. Conclusions

This paper presents nonlinear and linear methods for identifying the complex wavenumber space and estimating the DLF of heterogeneous metastructures in omni-direction. The methods are based on correlating the displacement field obtained by FEM with an inhomogeneous wave model (GFC) and an algebraic equation (AWI). The accuracy and robustness of the methods are numerically demonstrated by comparison with reference methods such as the WFE scheme and the PIM-FEM approach.

The WIM and the reference solution show good agreement in the dispersion curve and the DLF. Nonlinear and linear techniques can accurately predict the real part of the wavenumber over the entire frequency range. Both methods are also more reliable in HDS due to the free-field assumption where only the direct field is considered for wave identification.

In future research endeavors, analogous to the improvement of IWC, the accuracy of linear AWI can be theoretically optimized by employing Green’s function-based model instead of the exponential decay model, such development will theoretically optimize the accuracy of linear AWI.

Acknowledgements

The authors thank Dr. ZHOU Wei from Shenzhen University for his technical support in wave filter by physical constraints, and Dr. GUO Yunpeng for his technical support in conventional Finite Element Analysis.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Cui, D., Ichchou, M., Zine, A. and Atalla, N. (2024) Multi-Scale Dynamics and Nonlinear Eigenvalue Problem of Heterogeneous Metastructures Using a Wave Finite Element Scheme and Modal Strain Energy Method. In: Belhaq, M., Ed., Advances in Nonlinear Dynamics and Control of Mechanical and Physical Systems, Springer, 141-152.[CrossRef]
[2]	Cui, D., Atalla, N., Ichchou, M. and Zine, A. (2024) Damping Prediction of Highly Dissipative Meta-Structures through a Wave Finite Element Methodology. Mechanical Systems and Signal Processing, 215, Article ID: 111408.[CrossRef]
[3]	Cui, D., Ichchou, M., Atalla, N. and Zine, A. (2024) Computation of the Damping Loss Factor of Heterogeneous Meta-Structure Using the Wave Finite Element-Based Methodology. Journal of Physics: Conference Series, 2909, Article ID: 012008.[CrossRef]
[4]	Cui, D., Ichchou, M., Atalla, N. and Zine, A. (2025) Wave-Based Approaches for Wavespace of Highly Contrasted Structures with Viscoelastic Damping. Chinese Journal of Aeronautics, 38, Article ID: 103347.[CrossRef]
[5]	Cui, D., Ichchou, M., Atalla, N. and Zine, A. (2025) Physical Investigation on the Sound Transmission Loss of Heterogeneous Metastructures Using Wave-Based Methodologies. Aerospace Science and Technology, 161, Article ID: 110110.[CrossRef]
[6]	Lu, X. and Liu, D. (1992) Interlayer Shear Slip Theory for Cross-Ply Laminates with Nonrigid Interfaces. AIAA Journal, 30, 1063-1073.[CrossRef]
[7]	Sun, C. and Whitney, J.M. (1973) Theories for the Dynamic Response of Laminated Plates. AIAA Journal, 11, 178-183.[CrossRef]
[8]	Ford, R.D., Lord, P. and Walker, A.W. (1967) Sound Transmission through Sandwich Constructions. Journal of Sound and Vibration, 5, 9-21.[CrossRef]
[9]	Smolenski, C.P. and Krokosky, E.M. (1973) Dilational-Mode Sound Transmission in Sandwich Panels. The Journal of the Acoustical Society of America, 54, 1449-1457.[CrossRef]
[10]	Moore, J.A. and Lyon, R.H. (1991) Sound Transmission Loss Characteristics of Sandwich Panel Constructions. The Journal of the Acoustical Society of America, 89, 777-791.[CrossRef]
[11]	Narayanan, S. and Shanbhag, R.L. (1982) Sound Transmission through a Damped Sandwich Panel. Journal of Sound and Vibration, 80, 315-327.[CrossRef]
[12]	Ghinet, S. and Atalla, N. (2011) Modeling Thick Composite Laminate and Sandwich Structures with Linear Viscoelastic Damping. Computers & Structures, 89, 1547-1561.[CrossRef]
[13]	Zarastvand, M., Ghassabi, M. and Talebitooti, R. (2021) Prediction of Acoustic Wave Transmission Features of the Multilayered Plate Constructions: A Review. Journal of Sandwich Structures & Materials, 24, 218-293.[CrossRef]
[14]	Ghafouri, M., Ghassabi, M., Zarastvand, M.R. and Talebitooti, R. (2022) Sound Propagation of Three-Dimensional Sandwich Panels: Influence of Three-Dimensional Re-Entrant Auxetic Core. AIAA Journal, 60, 6374-6384.[CrossRef]
[15]	Atalla, N. and Sgard, F. (2015) Finite Element and Boundary Methods in Structural Acoustics and Vibration. Taylor & Francis Group.
[16]	Lodygowski, T. and Sumelka, W. (2006) Limitations in Application of Finite Element Method in Acoustic Numerical Simulation. Journal of Theoretical and Applied Mechanics, 44, 849-865. http://www.ptmts.org.pl/jtam/index.php/jtam/article/download/v44n4p849/497
[17]	Allard, J.F. and Atalla, N. (2009) Propagation of Sound in Porous Media: Modelling Sound Absorbing Materials. 2nd Edition, Wiley.[CrossRef]
[18]	Gupta, P. and Parey, A. (2022) Prediction of Sound Transmission Loss of Cylindrical Acoustic Enclosure Using Statistical Energy Analysis and Its Experimental Validation. The Journal of the Acoustical Society of America, 151, 544-560.[CrossRef] [PubMed]
[19]	Hambric, S.A., Hwang, Y.F. and Bonness, W.K. (2004) Vibrations of Plates with Clamped and Free Edges Excited by Low-Speed Turbulent Boundary Layer Flow. Journal of Fluids and Structures, 19, 93-110.[CrossRef]
[20]	Atalla, N. (2014) Modeling the Sound Transmission through Complex Structures with Attached Noise Control Materials. Wave Motion, 51, 650-663.[CrossRef]
[21]	De Rosa, S. and Franco, F. (2008) A Scaling Procedure for the Response of an Isolated System with High Modal Overlap Factor. Mechanical Systems and Signal Processing, 22, 1549-1565.[CrossRef]
[22]	Ichchou, M.N., Hiverniau, B. and Troclet, B. (2009) Equivalent ‘Rain on the Roof’ Loads for Random Spatially Correlated Excitations in the Mid-High Frequency Range. Journal of Sound and Vibration, 322, 926-940.[CrossRef]
[23]	Courtois, T., Bertolini, C. and Ochs, J. (2010) A Procedure for Efficient Trimmed Body FE Simulations, Based on a Transfer Admittance Model of the Sound Package. SAE International Journal of Passenger Cars—Mechanical Systems, 3, 1-13.[CrossRef]
[24]	Bin Fazail, M.N., Chazot, J., Lefebvre, G. and Atalla, N. (2023) Damping Loss Factor Characterization of Complex Structures Using a Green’s Function-Based Model. Journal of Sound and Vibration, 552, Article ID: 117642.[CrossRef]
[25]	Bin Fazail, M., Chazotand, J.D., Lefebvre, G. and Atalla, N. (2022) Identification of the Dispersion Curves and the Damping Loss Factor using Green’s Function-Based Model of Non-Isotropic Structures. Proceedings of the 28th International Congress on Sound and Vibration, Singapore, 24-28 July 2022, 1-8.
[26]	Li, X., Ichchou, M., Zine, A., Droz, C. and Bouhaddi, N. (2023) An Algebraic Wavenumber Identification (AWI) Technique under Stochastic Conditions. Mechanical Systems and Signal Processing, 188, Article ID: 109983.[CrossRef]
[27]	Tufano, G. (2020) K-Space Analysis of Complex Large-Scale Periodic Structures. Master’s Thesis, école Centrale de Lyon.
[28]	Papagiannopoulos, G.A. and Hatzigeorgiou, G.D. (2011) On the Use of the Half-Power Bandwidth Method to Estimate Damping in Building Structures. Soil Dynamics and Earthquake Engineering, 31, 1075-1079.[CrossRef]
[29]	Liu, C.Q. and Goetchius, G.M. (2001) Estimation of Damping Loss Factors by Using the Hilbert Transform and Exponential Average Method. SAE Technical Paper Series, 1, 1075-1079.[CrossRef]
[30]	Jaber, M., Schneeweiβ, H., Bös, J.B. and Melz, T. (2014) Measurement of the Damping Properties of Carbon Composite Plates by the Power Input Method. ISMA 2014-International Conference on Noise and Vibration Engineering, Leuven, 15-17 September 2014, 1445-1458.
[31]	Brion, T. (2024) Suivi de santé structurale sur structures architecturées par techniques en ondes guidées. Master’s Thesis, école Centrale de Lyon.
[32]	Berthaut, J. (2004) Contribution à l’identification large bande des structures anisotropes: Application aux tables d’harmonie des pianos. Master’s Thesis, école Centrale de Lyon.
[33]	Berthaut, J., Ichchou, M.N. and Jezequel, L. (2005) Space Identification of Apparent Structural Behaviour. Journal of Sound and Vibration, 280, 1125-1131.[CrossRef]
[34]	Cuenca, J., Gautier, F. and Simon, L. (2009) The Image Source Method for Calculating the Vibrations of Simply Supported Convex Polygonal Plates. Journal of Sound and Vibration, 322, 1048-1069.[CrossRef]
[35]	Gunda, R., Vijayakar, S.M., Singh, R. and Farstad, J.E. (1998) Harmonic Green’s Functions of a Semi-Infinite Plate with Clamped or Free Edges. The Journal of the Acoustical Society of America, 103, 888-899.[CrossRef]
[36]	Cremer, L., Heckl, M. and Ungar, E.E. (1988) Structure-Borne Sound: Structural Vibrations and Sound Radiation at Audio Frequencies. 2nd Edition, Springer.
[37]	Himmich, I., El Hallaoui, I. and Soumis, F. (2024) A Multiphase Dynamic Programming Algorithm for the Shortest Path Problem with Resource Constraints: A Guide to Numerical Dispersion Curve Calculations. European Journal of Operational Research, 315, 470-483.[CrossRef]
[38]	Alimonti, L., Atalla, N. and Gardner, B. (2018) Assessment of a Contour Integral Method for the Nonlinear Eigenvalue Problem Arising from the Wave Finite Element Description of Two-Dimensional Periodic Waveguides Involving Dissipative and Frequency Dependent Properties Media. The 6th Conference of Noise and Vibration Emerging Methods, Ibiza, 7-9 May 2018.
[39]	Reissner, E. (1945) The Effect of Transverse Shear Deformation on the Bending of Elastic Plates. Journal of Applied Mechanics, 12, A69-A77.[CrossRef]
[40]	Altenbach, H., Eremeyev, V.A. and Naumenko, K. (2015) On the Use of the First Order Shear Deformation Plate Theory for the Analysis of Three‐Layer Plates with Thin Soft Core Layer. ZAMM—Journal of Applied Mathematics and Mechanics/Zeitschrift für Angewandte Mathematik und Mechanik, 95, 1004-1011.[CrossRef]
[41]	Arasan, U., Marchetti, F., Chevillotte, F., Tanner, G., Chronopoulos, D. and Gourdon, E. (2021) On the Accuracy Limits of Plate Theories for Vibro-Acoustic Predictions. Journal of Sound and Vibration, 493, Article ID: 115848.[CrossRef]
[42]	Butaud, P., Foltête, E. and Ouisse, M. (2016) Sandwich Structures with Tunable Damping Properties: On the Use of Shape Memory Polymer as Viscoelastic Core. Composite Structures, 153, 401-408.[CrossRef]
[43]	Boutin, C., Viverge, K. and Hans, S. (2021) Dynamics of Contrasted Stratified Elastic and Viscoelastic Plates—Application to Laminated Glass. Composites Part B: Engineering, 212, Article ID: 108551.[CrossRef]