Local Peaks Search Method for Solving Lamb Waves’ Dispersion Equation of Laminated Structures and the Application

To study the acoustic characteristics of sound scattered from laminated structures such as elastic plates and shells, it is usually required to solve the Lamb waves’ dispersion equations. Many traditional root-finding methods such as bisection, the Newton–Raphson method, and the Muller method are not able to tackle the problem completely. A simple but powerful method named local peaks search (LPS) is proposed to overcome their drawbacks. Firstly, the non-zero part of the dispersion equation is defined as the dispersion function, and its reciprocal is used to transform the zeros (i.e., roots) into local peaks. Secondly, the chosen complex domain is discretized, and the coarse local domains where the local peaks exist are determined by the direct search method globally. Thirdly, the Muller method is applied to obtain the refined locations of local peaks. Lastly, in order to refine the results, a hierarchical scheme is designed and the iteration of the above procedures is implemented; the error is set to be 10−16 as the stop criteria. The accuracy of the LPS method is validated by comparing it with the bisection method for the problem of elastic plates in the vacuum. The acoustic echo structures are analyzed experimentally. By computation of Lamb waves’ phase velocity, the critical angles are derived numerically and compared with the results acquired by an experiment using monostatic sound transducers. In this way, it is validated that the elastic scattered wave components are the highlights shown in the time-angle figure. Furthermore, the work can be applied for non-destructive testing, especially underwater structural health monitoring.


Introduction
The acoustic characteristics of sound propagated in or scattered from plates, shells, and other laminated structures have strong relationships with Lamb waves.In underwater acoustics, Lamb waves are important to study the highlights of elastic echoes scattered from structures [1][2][3][4].Solving the dispersion equations is the foundational requirement; otherwise, the studies on acoustic characteristics are not able even to begin.
The problem of Lamb waves in a thin elastic plate with a liquid load was first studied by Osborne and Hart in 1945; the results could be expressed in terms of symmetric and anti-symmetric modes [1].Zhu and Wu considered the viscosity of the fluid and derived the dispersion equation for Lamb waves in thin plates with viscous fluid layers [2].Li et al. used the micro-perturbation method to approximate the load layer as appropriate boundary conditions and obtained a first-order approximate dispersion equation suitable for various thin film loads [3].Nayfeh and Nagy improved the work of Zhu and Wu by discussing three different models for viscous fluids.They calculated the dispersion equations for leaky Lamb waves in thin plates with viscous fluid layers under several different circumstances [4].Simonetti and Cawley studied the dispersion characteristics of SH waves in thin plates with viscoelastic layers [5,6].Lu et al. studied the propagation characteristics of Lamb waves in a thin plate with a viscoelastic thin layer on one side and proposed a method for determining the complex Lamé constants of viscoelastic films, which involves determining µ re by SH waves and determining µ im by Lamb waves [7].In these studies, there is widespread exploration of Lamb waves' applications; however, the method for solving the Lamb wave equations has not been explicitly elucidated, posing a significant challenge in deriving dispersion and attenuation curves when commencing research in this domain.
The dispersion equation of Lamb waves is a transcendental equation that cannot be solved analytically.Root-finding methods, such as the bisection method and Newton-Raphson method, are widely implemented to obtain the real roots [8,9].However, for many problems, the dispersion equations are complex, whose roots are also complex, which makes the problem much more difficult to solve.Barshinger and Rose proposed a minimization algorithm with real roots as the initial values to obtain the complex roots of the dispersion equation [10].Ross applied the iterative method to solve the real roots of the dispersion equation and introduced the Muller method for solving complex roots.Muller's method is proposed as a numerical root-finding algorithm for solving equations of the form f (z) = 0, where z can be real or complex.It was first presented by D. E. Muller in 1956 and is faster than the bisection method and overcomes one of the biggest drawbacks of the Newton-Raphson method, i.e., computation of derivatives at each step [11].Liu et al. used the function 'fminbnd' in MATLAB to find the minimum points based on the matrix representation of the dispersion equation [12].Some researchers used the winding number integral method to solve the complex roots of the dispersion equation.This method is based on the generalized residue theorem, which divides the complex plane into rectangular grids, determines the number of complex roots in the rectangular boundary first, and then combines algebraic theorems to obtain the complex roots of the equation [13,14].However, the winding number integral method is rather difficult to implement and not stable.By investigation, it can be found that these methods applied for analytical formulation of the Lamb waves' dispersion equation have their own drawbacks and are not be able to tackle with the problem completely, especially without missing any roots.
Recently, purely numerical methods such as the finite difference method, the finite element method, and the boundary element method have been studied by researchers.Uhl et al. used the finite difference method, which is a generalized semi-analytical method, to compute the dispersive curves.The method was proved stable for the considered frequency range (0~8 MHz•mm) [15].Lamb waves within a sandwich skin were simulated based on the finite element method, and ABAQUS Software was used to extract the dispersion curve [16].Semi-analytical formulations of Lamb waves and the dispersion curves were obtained by the scaled boundary finite element method (SBFEM), which combines advantages of the finite element and boundary element methods simultaneously [17].Pablo Gómez et al. used COMSOL Multiphysics to generate pseudo-experimental data and then used the MASW (multichannel analysis of surface waves) method to obtain the dispersion curves of Lamb waves in plates [18].However, these methods are difficult to implement and unstable and require a large number of elements, which makes the algorithm inefficient.Actually, these purely numerical methods are much more suitable for problems in which the analytic formulations of the dispersion equations cannot be derived.For Lamb waves in laminated structures, the analytical formulation of the dispersion equation can be expressed as a determinant.Therefore, root-finding methods are preferred for problems of Lamb waves concerning layered plates or laminated structures.However, since existing root-finding methods cannot deal with the problems perfectly, the LPS method is proposed to overcome their drawbacks while keeping their advantages.
In this paper, elastic plates are studied as the laminated structures.The local peaks search method is established, which does not require any characteristics of the dispersion function such as the derivatives, orthogonality, etc.The only requirement is that the absolute value of the complex function is continuous and can be computed accurately, which is the natural characteristic of the problem.In Section 2, the theory of Lamb waves is expressed in detail, and the dispersion equations are derived.In Section 3, the method of local peaks search is proposed, and the procedures are discussed thoroughly.In Section 4, the software is described in detail, and the algorithm of the method is further discussed.In Section 5, the LPS method is validated, and a problem of Lamb waves is computed; dispersion and attenuation curves are obtained.In Section 6, the structures of acoustic echoes scattered from an elastic plate are analyzed experimentally by monostatic sonars, and the critical angle is computed by the LPS method and compared with experimental results.

Theory
The problem of acoustic propagation in elastic bodies can be divided into two categories: bulk waves and guided waves, both of which are based on elastic theory and satisfy the same wave equation.Bulk waves are waves propagating in an infinitely elastic medium, without the constraints of boundary conditions, and can be decomposed into L-waves, SH-waves, and SV-waves.Guided waves are waves that propagate along a waveguide and are synthesized by bulk waves passing through boundary conditions.Therefore, certain boundary conditions must be met.It is the limitation of boundary conditions that makes the solution of guided wave problems very difficult.Some guided wave problems with simpler boundary conditions have been solved, such as Rayleigh waves, Lamb waves, Stoneley waves, Love waves, and Scholte waves.These guided waves are also known as surface waves or interface waves.
The problem shown in Figure 1 can be used to express many diverse boundary conditions for the Lamb waves in the elastic plate.The two boundary conditions of BC1 and BC2 can be any one of air, vacuum, and water and the like fluid materials.On the elastic plate, there usually exists a viscoelastic layer on one side, which is the usage of acoustic stealth.All of these different fluids or materials will lead to a difficult problem for solving the complex dispersion equation because the roots are all complex.Especially for the problem under diverse boundary conditions, a universal method is still needed.
In this paper, elastic plates are studied as the laminated structures.The local peaks search method is established, which does not require any characteristics of the dispersion function such as the derivatives, orthogonality, etc.The only requirement is that the absolute value of the complex function is continuous and can be computed accurately, which is the natural characteristic of the problem.In Section 2, the theory of Lamb waves is expressed in detail, and the dispersion equations are derived.In Section 3, the method of local peaks search is proposed, and the procedures are discussed thoroughly.In Section 4, the software is described in detail, and the algorithm of the method is further discussed.In Section 5, the LPS method is validated, and a problem of Lamb waves is computed; dispersion and attenuation curves are obtained.In Section 6, the structures of acoustic echoes scattered from an elastic plate are analyzed experimentally by monostatic sonars, and the critical angle is computed by the LPS method and compared with experimental results.

Theory
The problem of acoustic propagation in elastic bodies can be divided into two categories: bulk waves and guided waves, both of which are based on elastic theory and satisfy the same wave equation.Bulk waves are waves propagating in an infinitely elastic medium, without the constraints of boundary conditions, and can be decomposed into Lwaves, SH-waves, and SV-waves.Guided waves are waves that propagate along a waveguide and are synthesized by bulk waves passing through boundary conditions.Therefore, certain boundary conditions must be met.It is the limitation of boundary conditions that makes the solution of guided wave problems very difficult.Some guided wave problems with simpler boundary conditions have been solved, such as Rayleigh waves, Lamb waves, Stoneley waves, Love waves, and Scholte waves.These guided waves are also known as surface waves or interface waves.
The problem shown in Figure 1 can be used to express many diverse boundary conditions for the Lamb waves in the elastic plate.The two boundary conditions of BC1 and BC2 can be any one of air, vacuum, and water and the like fluid materials.On the elastic plate, there usually exists a viscoelastic layer on one side, which is the usage of acoustic stealth.All of these different fluids or materials will lead to a difficult problem for solving the complex dispersion equation because the roots are all complex.Especially for the problem under diverse boundary conditions, a universal method is still needed.For guided wave propagation in elastic plates, the displacement wave equation is used, which can be written as ( ) For guided wave propagation in elastic plates, the displacement wave equation is used, which can be written as in which, i, j, k = 1, 2, 3, which represents the x, y, and z axis.λ and µ are the Lamé constants.This equation can be used for elastic and viscoelastic materials.
In order to solve the problem, the potential function is used to derive the displacement, which has the following form: in which it can be seen that the displacement vector u can be decomposed into two components: one is the scaler potential Φ and the other is the vector potential Ψ.
For the problems shown in Figure 1, the strain is considered to be in a plane, thus the vector potential can be reduced to one dimension, that is, Ψ = Ψ 2 (x 1 , x 3 )j.In this way, the problem is largely simplified.The displacement potential function can be derived as in which it can be seen that the second term u 2 is zero.
Based on the knowledge of the plane problem, the strain can be written as In general, the solution of the equation can be represented by formulation of plane waves in which, Next, the formula solution of displacement is +βA t sin(βx 3 ) − βB t cos(βx 3 ) u 3 = −αA l sin(αx 3 ) + αB l cos(αx 3 ) +iξ A t cos(βx 3 ) + iξB t sin(βx 3 ) and the plane wave formula solution of the strain is in which the time factor exp(−i(ωt − ξx 1 )) is omitted.For viscoelastic materials, the complex Young's modulus is defined as ), E is the real Young's modulus, and η E is the loss factor; the relationship with complex Lame constants λ * and µ * can be derived, that is, The acoustic parameters of the viscoelastic media can be defined as the complex longitudinal wave speed and transverse wave speed, It can also introduce the attenuation coefficients as the acoustic parameters; that is, the complex longitudinal wavenumber and the complex transverse wavenumber are defined as Sensors 2023, 23, 9359 in which c L is the longitudinal wave speed, α L is the longitudinal attenuation coefficient, c T is the transverse wave speed, and α T is the transverse attenuation coefficient.Obviously, α L and α T are the image part of the corresponding wavenumber.
For the viscoelastic layer, the formal solution of the displacement potential function can be expressed as Substituted it into the wave equation, we can obtain In fluid, the wave equation of displacement potential function is under which u = ∇Φ can be assumed.The displacement and sound pressure are written as In fluid 1, assume the form solutions of displacement function, the displace, and the sound pressure are In fluid 2, the formal solutions of displacement function, the displace, and the sound pressure are expressed as Equations ( 6), (7), ( 12), ( 13), ( 16) and ( 17) are combined to derive the Lamb waves' dispersion equation under diverse boundary conditions, especially the very complicated one.Eventually, the dispersion equation can be written as where f (z) is denominated as the dispersion function.
To solve dispersion Equation (18), it is necessary to find the zeros of f (z).When the boundary condition becomes complicated, the dispersion equation will be rather difficult to solve.

Method
Since the bisection method and Muller method cannot calculate the complex roots of the complex dispersion equations with complicated diverse boundary conditions, a novel method nominated as the local peaks search (LPS) method is proposed, which overcomes the drawbacks of Muller method while keeping its advantages.The LPS method is also a root-finding method; it searches all the local peaks and utilizes the Muller method to obtain the solutions of the dispersion equations by using a hierarchical scheme.
The zeros of complex function f (z) can be solved by the LPS method by introducing a reciprocal function g(z): in which it can be seen that the zeros of f (z) are the poles of g(z).In Equation ( 19), the absolute value of f (z) is added by one to prevent overflow; that is, make sure the denominator is not zero.Consider the function g(z), which is analytic over the whole complex plane, excluding some poles.The function has the characteristic that its absolute value is infinite at the poles and decreases rapidly when z is away from the poles.Therefore, as shown in Figure 2, if drawing a 3D image over the complex plane, a local peak appears around the poles in the image.Through searching the locations of local peaks, the complex domain of poles can be determined coarsely.Furthermore, the accurate value of poles can be obtained by the iteration method implemented among the determined coarse domains.
Sensors 2023, 23, x FOR PEER REVIEW 7 of 23 Consider the function g(z), which is analytic over the whole complex plane, excluding some poles.The function has the characteristic that its absolute value is infinite at the poles and decreases rapidly when z is away from the poles.Therefore, as shown in Figure 2, if drawing a 3D image over the complex plane, a local peak appears around the poles in the image.Through searching the locations of local peaks, the complex domain of poles can be determined coarsely.Furthermore, the accurate value of poles can be obtained by the iteration method implemented among the determined coarse domains.Firstly, the global domain search procedure is implemented.The location of local peaks of g(z) can be determined by a direct search method over the global discretized zone.First, it should select the zone to search, that is, determine the size of the complex plane to be discretized empirically, which should not be too large or too small, and the step used for discretizing the real and imaginary parts should be sufficiently small, in order to guarantee that all local peaks are included.When searched all over the selected zone globally, the coarse rectangle domains in the complex plane where local peaks locate are found by the direct search method totally.
To find the local domains, it is needed to make sure there are local peaks inside.Assume the discretized values of g(z) are ( ) where the terms of window function ( ) Firstly, the global domain search procedure is implemented.The location of local peaks of g(z) can be determined by a direct search method over the global discretized zone.First, it should select the zone to search, that is, determine the size of the complex plane to be discretized empirically, which should not be too large or too small, and the step used for discretizing the real and imaginary parts should be sufficiently small, in order to guarantee that all local peaks are included.When searched all over the selected zone globally, the coarse rectangle domains in the complex plane where local peaks locate are found by the direct search method totally.
To find the local domains, it is needed to make sure there are local peaks inside.Assume the discretized values of g where the terms of window function w(i, j) are shown in Figure 3. sume the discretized values of g(z) are ( ) where the terms of window function ( ) , w i j are shown in Figure 3.
If g(m, n) is the maximum of the 9 values inside the window, then g(m, n) might be the peak value of a local peak.When far away from the zeros, the value of the function changes slowly, while near the zeros, the value of the function changes dramatically.Using this characteristic, it is possible to determine whether it is a local peak or a small bump based on the speed of change of the function value.A metric ∆(m, n) is proposed, which can be used to measure the change speed of function value: The reason for adopting this metric is that the value of g(z) varies rather rapidly close to the poles of g(z); however, if far away from the poles, the value will change smoothly.Based on this characteristic, ∆(m, n) can judge the slope of the function value changes.A parameter of ε > 0 is chosen empirically; when ∆(m, n) ≥ ε, the local peaks could be determined to exist, as seen in Figure 2. In the following, the kurtosis of ε is set to be 0.001.This value can be chosen empirically or experimentally.
When local domains containing local peaks are found, the local domain search procedure is implemented.The local peaks found are not very accurate, since they are just the approximate of the true results.In order to promote accuracy, the local peak is used as the centroid to establish a rectangle domain whose size is about several steps for discretizing the real and imaginary parts of the complex plane, as shown in Figure 4. ( ) , The reason for adopting this metric is that the value of g(z) varies rather rapidly close to the poles of g(z); however, if far away from the poles, the value will change smoothly.Based on this characteristic, Δ(m, n) can judge the slope of the function value changes.A parameter of ε > 0 is chosen empirically; when ( ) ≥ , the local peaks could be determined to exist, as seen in Figure 2. In the following, the kurtosis of ε is set to be 0.001.This value can be chosen empirically or experimentally.When local domains containing local peaks are found, the local domain search procedure is implemented.The local peaks found are not very accurate, since they are just the approximate of the true results.In order to promote accuracy, the local peak is used as the centroid to establish a rectangle domain whose size is about several steps for discretizing the real and imaginary parts of the complex plane, as shown in Figure 4. Next, to find the accurate solution of the dispersion equation, the Muller method is used to compute the value of the roots within the rectangle domain.The Muller method is suitable to solve the equation with one real or complex root inside a contour (i.e., the rectangle domain in Figure 4).The reason to use the Muller method is that, while keeping the least requirement for the dispersion function, its precision can be rather high.Thus, the discretization step has to be small enough to make every grid contain only one local peak.
Finally, in order to guarantee the accuracy, a hierarchical scheme is established in the local domains; the iteration method is used to execute the scheme.When the coarse root is solved by Muller method, a refined rectangle domain around the root is set.Then the domain is discretized and the Muller method is used to refine the result again and again to make the error to be less than a given one, e.g., 10 −16 , which is also the stopping criteria of the iteration.In this way, a highly accurate result of the root can be obtained eventually.Next, to find the accurate solution of the dispersion equation, the Muller method is used to compute the value of the roots within the rectangle domain.The Muller method is suitable to solve the equation with one real or complex root inside a contour (i.e., the rectangle domain in Figure 4).The reason to use the Muller method is that, while keeping the least requirement for the dispersion function, its precision can be rather high.Thus, the discretization step has to be small enough to make every grid contain only one local peak.
Finally, in order to guarantee the accuracy, a hierarchical scheme is established in the local domains; the iteration method is used to execute the scheme.When the coarse root is solved by Muller method, a refined rectangle domain around the root is set.Then the domain is discretized and the Muller method is used to refine the result again and again to make the error to be less than a given one, e.g., 10 −16 , which is also the stopping criteria of the iteration.In this way, a highly accurate result of the root can be obtained eventually.

Software
The source code of this work is fully tested on various laminated structures under many boundary conditions, which are available on Github (see 'Data Availability Statement').In this section, the algorithm and the software are thoroughly described for usage.The source code is written in MATLAB programming language and tested on the R2022b version.
Figure 5 is the 'main' page of the software; it sets the parameters for the LPS method.For the 'Problem Selection' module, there are two different laminated structures to choose: one is the planar plate problem, and another is the cylindrical shell problem.Even though, in this paper, the planar plate problem is studied, the LPS method can also be applied for the cylindrical shell problems.The module of 'Mode Selection' is used only for cylindrical shell problems.For the 'Model Size' module, the thickness of the elastic and viscoelastic plates can be set; a, b, and c are the thickness of the elastic or viscoelastic cylindrical shells.For the 'Material Parameters' module, the parameters of fluid 1 and 2, elastic and viscoelastic materials are set.And the last module is 'Computation Parameters', in which the frequency range, the start, stop, and step parameters of the phase velocity, and the attenuation coefficient are set for discretization range of the complex plane.The 'Control Parameter' is set for determining the local peaks and the solution precision.The architecture of the algorithm is shown in Figure 6; the algorithm is an implementation of the method depicted in Section 3. Through the main function, all parameters needed are initialized.Then, the global domain search procedure is executed, which first discretizes the selected domain over the complex plane with the phase velocity and the attenuation coefficient as the x and y axis, respectively.The result of discretization can be seen in Figure 2, which contains only one local peak.To find all of the local peaks, the direct search method is used, which is called 'LocalPeaksSearch' in the algorithm.After that, the roots-find method is applied here to obtain a relatively more accurate result for further usage.In this way, all of the local peaks in the region of interest are found.The architecture of the algorithm is shown in Figure 6; the algorithm is an implementation of the method depicted in Section 3. Through the main function, all parameters needed are initialized.Then, the global domain search procedure is executed, which first discretizes the selected domain over the complex plane with the phase velocity and the attenuation coefficient as the x and y axis, respectively.The result of discretization can be seen in Figure 2, which contains only one local peak.To find all of the local peaks, the direct search method is used, which is called 'LocalPeaksSearch' in the algorithm.After that, the roots-find method is applied here to obtain a relatively more accurate result for further usage.In this way, all of the local peaks in the region of interest are found.Next, the local domain search method is executed; the procedures are the same as the global domain search.However, three procedures, i.e., the dispersion equation discretization, local peaks search, and the roots finding, are applied on the local domains around the local peaks.These three procedures are implemented to every single local peak.The error is controlled within the module 'RootsFind'.
In order to obtain a highly accurate solution, a mixed scheme is designed.Along with the root-find method based on the Muller method, a domain refine method is proposed based on the hierarchical scheme in Section 3.That is, as we continue to discretize the local domain containing local peaks, its size will be reduced iteratively to less than the error.In this way, the roots can also be obtained accurately.The mixed scheme is then designed to  Next, the local domain search method is executed; the procedures are the same as the global domain search.However, three procedures, i.e., the dispersion equation discretization, local peaks search, and the roots finding, are applied on the local domains around the local peaks.These three procedures are implemented to every single local peak.The error is controlled within the module 'RootsFind'.
In order to obtain a highly accurate solution, a mixed scheme is designed.Along with the root-find method based on the Muller method, a domain refine method is proposed based on the hierarchical scheme in Section 3.That is, as we continue to discretize the local domain containing local peaks, its size will be reduced iteratively to less than the error.In this way, the roots can also be obtained accurately.The mixed scheme is then designed to combine the Muller root-find method and the domain refine method, which makes the algorithm much more robust.
After the local domain search is executed, if there are multiple roots, a scheme is established to remove all of the repeated ones.Then, a manual check module can be chosen to check whether the root is right or not.Of course, this manual procedure is not required and can be cancelled at the beginning.Finally, the post processing module can be used to plot the dispersion and attenuation curves.

Validation
For the computation and analysis of the Lamb waves' dispersion curves, several different boundary conditions are selected to validate the LPS method.The material of the elastic plate is steel, its density is × 10 3 kg/m 3 , Young's modulus is 3.0 × 10 7 Pa, and the Poisson ratio is 0.28.
As shown in Figure 7, for Lamb waves' dispersion equations of free elastic plates, the bisection method can be used to compute the roots.The LPS method is compared with the bisection method in Figure 7; it can be seen that the results of the two methods are the same.Therefore, the accuracy of the LPS method is validated; it is highly accurate.It is absolutely a challenge to find a method to solve the complex dispersion equation of Lamb waves under diverse complicated boundary conditions.The bisection method failed in this problem, the LPS method, is applied to compute the dispersion and attenuation curves.To prove the potential of the LPS method to tackle such difficulties robustly, the dispersion equation of elastic plates coated with a viscoelastic layer in the vacuum is solved, as shown in Figure 8.The boundary conditions of the problem shown in Figure 8 can be written as  By using the boundary conditions, the dispersion equation can be obtained: in which the terms of the dispersion equation are given in Appendix A.
In the problem of Figure 8, the dispersion equation of the elastic plate coated with a viscoelastic damping layer in the vacuum is computed by the LPS method.The results of dispersion curves and attenuation curves are shown in Figure 9.In this problem, the symmetric and anti-symmetric modes of Lamb waves cannot be split; the first eight modes are computed.Furthermore, even though the modes cannot be split, the dispersion curves are somewhat like Figure 7, which means the viscoelastic layer is a load to the elastic plate and mainly influences the energy attenuation.The reason is that the impedance of the steel plate is much bigger than that of the viscoelastic layer.

Highlight Model
For the active sonar, elastic scattering echoes are closely related to the propagation of acoustic waves in elastic bodies, so studying the propagation of acoustic waves in elastic plates and cylindrical shells can help in the understanding of the structure of target echoes, as shown in Figure 10.When a plane wave is incident on an elastic plate, a portion of the acoustic wave is reflected from the surface of the elastic plate to form a specular reflection wave.A part of the sound wave is scattered by the edges and corners of the elastic plate to form edge waves; another part of the sound wave penetrates into the elastic plate and excites the formation of Lamb waves.Due to the coupling between water and the By using the boundary conditions, the dispersion equation can be obtained: in which the terms of the dispersion equation are given in Appendix A.
In the problem of Figure 8, the dispersion equation of the elastic plate coated with a viscoelastic damping layer in the vacuum is computed by the LPS method.The results of dispersion curves and attenuation curves are shown in Figure 9.In this problem, the symmetric and anti-symmetric modes of Lamb waves cannot be split; the first eight modes are computed.Furthermore, even though the modes cannot be split, the dispersion curves are somewhat like Figure 7, which means the viscoelastic layer is a load to the elastic plate and mainly influences the energy attenuation.The reason is that the impedance of the steel plate is much bigger than that of the viscoelastic layer.

Application 6.1. Highlight Model
For the active sonar, elastic scattering echoes are closely related to the propagation of acoustic waves in elastic bodies, so studying the propagation of acoustic waves in elastic plates and cylindrical shells can help in the understanding of the structure of target echoes, as shown in Figure 10.When a plane wave is incident on an elastic plate, a portion of the acoustic wave is reflected from the surface of the elastic plate to form a specular reflection wave.A part of the sound wave is scattered by the edges and corners of the elastic plate to form edge waves; another part of the sound wave penetrates into the elastic plate and excites the formation of Lamb waves.Due to the coupling between water and the elastic plate, the energy of the Lamb waves propagating in the plate leaks into the water, forming elastic scattered echoes, which is called the highlight model [19].
Sensors 2023, 23, x FOR PEER REVIEW 13 of 23 elastic plate, the energy of the Lamb waves propagating in the plate leaks into the water, forming elastic scattered echoes, which is called the highlight model [19].Considering elastic scattering echoes, the phase velocity of Lamb waves propagating in an elastic plate is denoted as p c , while the phase velocity of acoustic waves in water is w c .The boundary condition at the interface between water and the elastic plate requires the wavenumber of the acoustic wave to satisfy the continuity condition.Assuming that the angle between the elastic scattering wave and the normal of the boundary condition is Along the direction of angle c θ , the received echo signal contains the elastic scatter- ing echoes, while no elastic scattering wave is received in adjacent directions, resulting in a significant enhancement of the echo signal in the direction of c θ , known as the highlight points.c θ is called the critical angle.Since Lamb waves are multi-modal, different modes have different critical angles for exiting into acoustic waves in water, leading to the observation of elastic highlights at multiple different angles.

The Experiment Deployment
The experiment was conducted in a tank, whose size is 2m 1.5m 1.5m × × . Figure 11 shows the layout of the experimental instruments and equipment used in the experiment, as well as the arrangement of the transducer and the target.The transducer was installed using a transceiver configuration, and the target was suspended from a mechanic rotating device at the same depth as the transducer using a thin rope.The experiment is conducted by rotating the mechanical rotating device each step increased by one degree.The transmitting transducer emits a CW pulse, which is scattered by the elastic or laminated plate.The scattered waves are received by the receiving transducer.The mechanical device is rotated from 0 to 360 degrees, and experimental results are plotted in the time-angle pseudo-color image from which the highlights can be seen and determined.Considering elastic scattering echoes, the phase velocity of Lamb waves propagating in an elastic plate is denoted as c p , while the phase velocity of acoustic waves in water is c w .The boundary condition at the interface between water and the elastic plate requires the wavenumber of the acoustic wave to satisfy the continuity condition.Assuming that the angle between the elastic scattering wave and the normal of the boundary condition is Along the direction of angle θ c , the received echo signal contains the elastic scattering echoes, while no elastic scattering wave is received in adjacent directions, resulting in a significant enhancement of the echo signal in the direction of θ c , known as the highlight points.θ c is called the critical angle.Since Lamb waves are multi-modal, different modes have different critical angles for exiting into acoustic waves in water, leading to the observation of elastic highlights at multiple different angles.

The Experiment Deployment
The experiment was conducted in a tank, whose size is 2 m × 1.5 m × 1.5 m. Figure 11 shows the layout of the experimental instruments and equipment used in the experiment, as well as the arrangement of the transducer and the target.The transducer was installed using a transceiver configuration, and the target was suspended from a mechanic rotating device at the same depth as the transducer using a thin rope.The experiment is conducted by rotating the mechanical rotating device each step increased by one degree.The transmitting transducer emits a CW pulse, which is scattered by the elastic or laminated plate.The scattered waves are received by the receiving transducer.The mechanical device is rotated from 0 to 360 degrees, and experimental results are plotted in the time-angle pseudo-color image from which the highlights can be seen and determined.
As shown in Figure 11d, a set of self-developed transducers in the laboratory, labeled A1 and A2, were used.One of these transducers was used as the transmitting transducer and the other as the receiving transducer.The receiving transducer was rotated to measure its directivity.The frequencies of the transmitted signals in the experiment were set at 800 kHz and 1.1 MHz, and the directivity obtained can be seen in Figure 12.At 800 kHz, the −3 dB beamwidth of transducer A1 was approximately 5.2 • , and at 1.1 MHz, the −3 dB beamwidth was approximately 4.4 • .For transducer A2, at 800 kHz the −3 dB beamwidth was approximately 8.2 • , and at 1.1 MHz it was approximately 4.4 • .From the directivity measurements of A1 and A2, it is evident that the directivity of A2 is superior to that of A1.Therefore, when conducting the target echo experiment, it is advisable to select A1 as the transmitter and A2 as the receiver.As shown in Figure 11d, a set of self-developed transducers in the laboratory, labeled A1 and A2, were used.One of these transducers was used as the transmitting transducer and the other as the receiving transducer.The receiving transducer was rotated to measure its directivity.The frequencies of the transmitted signals in the experiment were set at 800 kHz and 1.1 MHz, and the directivity obtained can be seen in Figure 12.At 800 kHz, the −3 dB beamwidth of transducer A1 was approximately 5.2  , and at 1.1 MHz, the −3 dB beamwidth was approximately 4.4  .For transducer A2, at 800 kHz the −3 dB beamwidth was approximately 8.2  , and at 1.1 MHz it was approximately 4.4  .From the directivity measurements of A1 and A2, it is evident that the directivity of A2 is superior to that of A1.Therefore, when conducting the target echo experiment, it is advisable to select A1 as the transmitter and A2 as the receiver.As shown in Figure 11d, a set of self-developed transducers in the laboratory, labeled A1 and A2, were used.One of these transducers was used as the transmitting transducer and the other as the receiving transducer.The receiving transducer was rotated to measure its directivity.The frequencies of the transmitted signals in the experiment were set at 800 kHz and 1.1 MHz, and the directivity obtained can be seen in Figure 12.At 800 kHz, the −3 dB beamwidth of transducer A1 was approximately 5.2  , and at 1.1 MHz, the −3 dB beamwidth was approximately 4.4  .For transducer A2, at 800 kHz the −3 dB beamwidth was approximately 8.2  , and at 1.1 MHz it was approximately 4.4  .From the directivity measurements of A1 and A2, it is evident that the directivity of A2 is superior to that of A1.Therefore, when conducting the target echo experiment, it is advisable to select A1 as the transmitter and A2 as the receiver.

Case for the Steel Plate
The first target is a steel plate, named as steel plate 1; its size was 150.02 131.82 1.07 × × m; its density was .Figure 13 is the problem of the Lamb waves' dispersion equation of the elastic plate in water; the boundary condition can be described as

Case for the Steel Plate
The first target is a steel plate, named as steel plate 1; its size was 150.02 × 131.82 × 1.07 m; its density was ρ = 7.84 × 10 3 kg/m 3 ; Young's modulus was E = 21.6 × 10 10 N/m 2 ; Poisson's coefficient was σ = 0.28. Figure 13 is the problem of the Lamb waves' dispersion equation of the elastic plate in water; the boundary condition can be described as Under this boundary condition, the dispersion equation can be derived with the formal solutions.It can be found that, in this problem, just the same as the problem of elastic plates in vacuum, the vibration can be split into two different modes, i.e., the symmetric and anti-symmetric modes.The dispersion equation of the symmetric mode is and the dispersion equation of the anti-symmetric mode is

Case for the Steel Plate
The first target is a steel plate, named as steel plate 1; its size was 150.02 131.82 1.07 × × m; its density was .Figure 13 is the problem of the Lamb waves' dispersion equation of the elastic plate in water; the boundary condition can be described as Under this boundary condition, the dispersion equation can be derived with the formal solutions.It can be found that, in this problem, just the same as the problem of elastic plates in vacuum, the vibration can be split into two different modes, i.e., the symmetric and anti-symmetric modes.The dispersion equation of the symmetric mode is and the dispersion equation of the anti-symmetric mode is The dispersion equation of the steel plate in water is computed, as shown in Figure 14.In this problem, the water is served as the load to the plates, which will make the Lamb The dispersion equation of the steel plate in water is computed, as shown in Figure 14.In this problem, the water is served as the load to the plates, which will make the Lamb waves leaky; an attenuation coefficient is introduced to express the energy loss when Lamb waves propagate along the plate.However, the problem is symmetric; the Lamb waves can be split into symmetric and anti-symmetric modes.In this case, the dispersion function is complex, which results in complex roots of the solution.Using the LPS method to solve the problem, A0~A4 are the first five anti-symmetric modes and S0~S4 represent the first five symmetric modes.waves leaky; an attenuation coefficient is introduced to express the energy loss when Lamb waves propagate along the plate.However, the problem is symmetric; the Lamb waves can be split into symmetric and anti-symmetric modes.In this case, the dispersion function is complex, which results in complex roots of the solution.Using the LPS method to solve the problem, A0~A4 are the first five anti-symmetric modes and S0~S4 represent the first five symmetric modes.
(a) (b) From Figure 15, it can be seen that there are several highlights in the time-angle figure, which is related to the Lamb waves by the critical angle.In Table 1, the numerical and experimental critical angles are compared, which demonstrates that the highlights are related to the elastic scattering wave components of the echoes.The elastic scattering wave components are coming from the leakage of Lamb waves, which is the mechanism of the From Figure 15, it can be seen that there are several highlights in the time-angle figure, which is related to the Lamb waves by the critical angle.In Table 1, the numerical and experimental critical angles are compared, which demonstrates that the highlights are related to the elastic scattering wave components of the echoes.The elastic scattering wave components are coming from the leakage of Lamb waves, which is the mechanism of the acoustic echoes.From Figure 15, it can be seen that there are several highlights in the time-angle figure, which is related to the Lamb waves by the critical angle.In Table 1, the numerical and experimental critical angles are compared, which demonstrates that the highlights are related to the elastic scattering wave components of the echoes.The elastic scattering wave components are coming from the leakage of Lamb waves, which is the mechanism of the acoustic echoes.In Table 1, the phase speed of the Lamb wave is computed in the case shown in Figure 14.Using the formula of Equation ( 24), the critical angle is computed.From Figure 15, the critical angle can be determined by the highlights experimentally.Through comparison of the critical angle between numerical and experimental results, it can be seen the two results are close.This demonstrates the theory of the highlight model.In Table 1, the phase speed of the Lamb wave is computed in the case shown in Figure 14.Using the formula of Equation ( 24), the critical angle is computed.From Figure 15, the critical angle can be determined by the highlights experimentally.Through comparison of the critical angle between numerical and experimental results, it can be seen the two results are close.This demonstrates the theory of the highlight model.Another experiment is also conducted on a steel plate, named steel plate 2, whose size is 151.10 mm × 133.20 mm × 1.87 mm.After measuring the echoes by rotating the mechanic device from 0 to 180 degrees increased by one degree each step, the time-angle figure is obtained and drawn in Figure 16.The highlights can also be found in Figure 16a,b, the critical angles are determined, as shown in Table 2.The numerical results are also listed in Table 2, which are computed by the LPS method. .After measuring the echoes by rotating the mechanic device from 0 to 180 degrees increased by one degree each step, the time-angle figure is obtained and drawn in Figure 16.The highlights can also be found in Figure 16a,b, the critical angles are determined, as shown in Table 2.The numerical results are also listed in Table 2, which are computed by the LPS method., and the thickness is 1.87 mm.For the viscoelastic material, its density is 1.1 × 10 3 kg/m 3 , Young's modulus is 21.6 × 10 10 Pa, loss factor is 0.23, and Poisson's ratio is 0.49; the thickness is 2.39 mm.

Case for the Steel Plate Coated with a Viscoelastic Layer
In Figure 17, the material of the elastic plate is also steel, and it is coated with a viscoelastic layer, whose material is polyurethane.The size of the steel plate is 150.02 mm × 131.82 mm × 4.26 mm, and the thickness is 1.87 mm.For the viscoelastic material, its density is 1.1 × 10 3 kg/m 3 , Young's modulus is 21.6 × 10 10 Pa, loss factor is 0.23, and Poisson's ratio is 0.49; the thickness is 2.39 mm.For the dispersion equation of the elastic coated with a viscoelastic damping layer under water, the boundary condition can be written as

h d x h d x h d x h d x h d x h x h x h x h x h x h x h x h x
By applying this boundary condition, the dispersion equation is derived: in which the terms of the equation are given in Appendix B. This case is more complicated; however, the water is still a load to the elastic plates.In this problem, the coated plate is immersed in water, which can be further used to analyze the acoustic characteristics of elastic echoes scattered from the structure of plates in the water, as shown in Figure 18.For the dispersion equation of the elastic coated with a viscoelastic damping layer under water, the boundary condition can be written as By applying this boundary condition, the dispersion equation is derived: in which the terms of the equation are given in Appendix B. This case is more complicated; however, the water is still a load to the elastic plates.In this problem, the coated plate is immersed in water, which can be further used to analyze the acoustic characteristics of elastic echoes scattered from the structure of plates in the water, as shown in Figure 18.This case is more complicated; however, the water is still a load to the elastic plates.In this problem, the coated plate is immersed in water, which can be further used to analyze the acoustic characteristics of elastic echoes scattered from the structure of plates in the water, as shown in Figure 18.From Figure 19, the critical angle can be determined by the highlights experimentally.The phase speed of the Lamb wave is computed in Figure 18, and then the critical angle is calculated afterward.In Table 3, through comparison of the critical angles between numerical and experimental results, it can further demonstrate the theory of the highlight model.However, some of the errors are large, caused by many factors that are analyzed in the next section.From Figure 19, the critical angle can be determined by the highlights experimentally.The phase speed of the Lamb wave is computed in Figure 18, and then the critical angle is calculated afterward.In Table 3, through comparison of the critical angles between numerical and experimental results, it can further demonstrate the theory of the highlight model.However, some of the errors are large, caused by many factors that are analyzed in the next section.

Error Analysis
Tables 1-3 present the theoretical calculations and experimental measurements of the critical angles for the highlights of elastic scattering echoes for steel plate 1, steel plate 2, and the coated steel plate.It can be observed from the corresponding results in the tables that the theoretical calculations and experimental measurements are generally consistent, with a certain degree of error in the angles between the two.However, the errors are evidently within the acceptable range for engineering purposes.The errors between numerical and experimental critical angles can be several degrees.Some reasons for the errors are as follows: 1.
The rotation device lacks precision, resulting in significant calibration errors on the turntable.When rotating one full circle, the actual rotation angle is approximately 9.73 • , whereas theoretically it should be 10 • .Additionally, the experiment only allows for rotation in increments of one large division (approximately 0.973 • degrees, with 10 divisions in one full rotation), which lacks precision in angle adjustment; 2.
Ensuring the accuracy of the corresponding angles in the experiment is challenging.For instance, the reflection point of the mirror should appear at angle 90 • , but there are deviations observed in the experiment.This necessitates a certain degree of angle correction during the data processing, introducing errors in determining the critical angle value; 3.
Due to the acoustic signal emitted by the transducer having a certain beamwidth, the elastic scattering highlights in the pseudo-color figure also exhibit some angular broadening, posing significant challenges in determining the critical angle size; 4.
The power amplifier used in the experiment has some issues, resulting in signal distortion and inadequate power output, causing the transducer to generate weak acoustic waves and insufficient signal-to-noise ratio; 5.
The material parameters used in the theoretical calculations are imprecise, leading to inaccurate calculation results that deviate from the measured experimental values.

Conclusions
This study on acoustic characteristics of sound propagated in or scattered from elastic plates and other laminated structures are related to Lamb waves.The key is to solve the dispersion equations of Lamb waves.Traditional root-finding methods are not able to tackle with the problem of solving dispersion equations of Lamb waves in the laminated structures totally, so a novel method named the LPS method is proposed to overcome their drawbacks while keeping their advantages.
The LPS method does not require computing the dispersion function's derivatives or its orthogonality, etc.; the only requirement is that the absolute value of the dispersion function is continuous, which is the natural characteristic for the problem.Therefore, the LPS method might be considered as a universal root-finding method to solve problems of laminated structures whose dispersion equations can be expressed analytically or even any real and complex equation f (z) = 0 if the absolute value of f (z) is continuous.
The work provides a fundamental tool for non-destructive testing, especially for underwater or industrial structural health monitoring.In this paper, the LPS software is applied to analyze the characteristics of the sound scattering from laminated structures of plates, i.e., the highlights formed by Lamb waves.The echoes are plotted in the timeangle figure to find and analyze the highlights and obtain the critic angle.Through computation of Lamb waves' phase velocity, the critical angle is computed and compared with the experimental results, which can be used to predict the direction of highlights.The comparison reveals the mechanism of the elastic scattering wave components of the acoustic echoes, which are formed by the Lamb waves leaking into water.The correctness of the mechanism for the formation of elastic scattering echoes explained by the Lamb theory was validated by comparing the results of theoretical calculations at the critical angle with experimental measurements.

Figure 1 .
Figure 1.The Lamb wave model for elastic plates under diverse complicated boundary conditions.

Figure 1 .
Figure 1.The Lamb wave model for elastic plates under diverse complicated boundary conditions.

Figure 2 .
Figure 2. A local peak of g(z) formed by the zero of f(z).

Figure 2 .
Figure 2. A local peak of g(z) formed by the zero of f (z).

Sensors 2023 ,
23, x FOR PEER REVIEW 8 of 23 this characteristic, it is possible to determine whether it is a local peak or a small bump based on the speed of change of the function value.A metric Δ(m, n) is proposed, which can be used to measure the change speed of function value:

Figure 4 .
Figure 4. Rectangle domain containing the local peak.

Figure 4 .
Figure 4. Rectangle domain containing the local peak.

Sensors 2023 ,
23, x FOR PEER REVIEW 9 of 23 one is the planar plate problem, and another is the cylindrical shell problem.Even though, in this paper, the planar plate problem is studied, the LPS method can also be applied for the cylindrical shell problems.The module of 'Mode Selection' is used only for cylindrical shell problems.For the 'Model Size' module, the thickness of the elastic and viscoelastic plates can be set; a, b, and c are the thickness of the elastic or viscoelastic cylindrical shells.For the 'Material Parameters' module, the parameters of fluid 1 and 2, elastic and viscoelastic materials are set.And the last module is 'Computation Parameters', in which the frequency range, the start, stop, and step parameters of the phase velocity, and the attenuation coefficient are set for discretization range of the complex plane.The 'Control Parameter' is set for determining the local peaks and the solution precision.

Figure 5 .
Figure 5.The main page for the LPS method based on MATLAB.

Figure 5 .
Figure 5.The main page for the LPS method based on MATLAB.

Figure 6 .
Figure 6.The algorithm of the LPS method program.

Figure 6 .
Figure 6.The algorithm of the LPS method program.

Corrigendum 1 、
Figure 7.Comparison of the LPS method and bisection method.

2 、Figure 7 .
Figure 7.Comparison of the LPS method and bisection method.

Figure 8 .
Figure 8. Lamb wave model of elastic plates coated with a viscoelastic layer in the vacuum.

Figure 9 .
Figure 9. Dispersion and attenuation curves of Lamb waves from elastic plates with a viscoelastic layer in the vacuum.(a) Dispersion curve; (b) attenuation curve.

Figure 8 .
Figure 8. Lamb wave model of elastic plates coated with a viscoelastic layer in the vacuum.

Corrigendum 1 、
The vertical axis label of Figure7is wrong, the corrected one is list below:

Figure 7 .
Figure 7.Comparison of the LPS method and bisection method.

2 、Figure 9 .
Figure 9. Dispersion and attenuation curves of Lamb waves from elastic plates with a viscoelastic layer in the vacuum.(a) Dispersion curve; (b) attenuation curve.

Figure 9 .
Figure 9. Dispersion and attenuation curves of Lamb waves from elastic plates with a viscoelastic layer in the vacuum.(a) Dispersion curve; (b) attenuation curve.

Figure 10 .
Figure 10.Echo structure of elastic plates in water.

Figure 10 .
Figure 10.Echo structure of elastic plates in water.

Figure 11 .
Figure 11.Diagram of connection and layout of experimental instruments and equipment.(a) The instruments for sound transmitting and receiving; (b) the mechanical equipment for deployment; (c) the deployment of the elastic plate; (d) the deployment of the monostatic sonar transducers (acoustic sensors).

Figure 11 .
Figure 11.Diagram of connection and layout of experimental instruments and equipment.(a) The instruments for sound transmitting and receiving; (b) the mechanical equipment for deployment; (c) the deployment of the elastic plate; (d) the deployment of the monostatic sonar transducers (acoustic sensors).

Figure 11 .
Figure 11.Diagram of connection and layout of experimental instruments and equipment.(a) The instruments for sound transmitting and receiving; (b) the mechanical equipment for deployment; (c) the deployment of the elastic plate; (d) the deployment of the monostatic sonar transducers (acoustic sensors).

Figure 13 .
Figure 13.Leaky Lamb wave mode of elastic plates in water.

Figure 13 .
Figure 13.Leaky Lamb wave mode of elastic plates in water.

Figure 14 .
Figure 14.Dispersion and attenuation curves of Lamb waves from elastic plates in water.(a) Dispersion curve; (b) attenuation curve.

Figure 14 .
Figure 14.Dispersion and attenuation curves of Lamb waves from elastic plates in water.(a) Dispersion curve; (b) attenuation curve.

Figure 14 .
Figure 14.Dispersion and attenuation curves of Lamb waves from elastic plates in water.(a) Dispersion curve; (b) attenuation curve.

Figure 15 .
Figure 15.Time-angle pseudo-color picture of acoustic echoes from the steel plate.(a) Frequency 800 kHz; (b) frequency 1.1 MHz.

Figure 15 .
Figure 15.Time-angle pseudo-color picture of acoustic echoes from the steel plate.(a) Frequency 800 kHz; (b) frequency 1.1 MHz.

Sensors 2023 ,Figure 17 .
Figure 17.Lamb wave model for the elastic plate coated with viscoelastic layer in fluid.

Fluid 1 :Figure 17 .
Figure 17.Lamb wave model for the elastic plate coated with viscoelastic layer in fluid.
terms of the equation are given in Appendix B.

Figure 18 .
Figure 18.Dispersion and attenuation curves of Lamb waves from elastic plates coated with viscoelastic layer in water.(a) Dispersion curve; (b) attenuation curve.As shown in Figure 19, the data of acoustic scattering echoes is measured from the side of steel; the time-angle figure is obtained by experiment for different frequencies of 800 kHz and 1.1 MHz.From Figures 15, 16 and 19, it is possible to clearly observe the highlights of mirror reflection echoes, edge wave echoes, and elastic scattering echoes.

Figure 18 .Figure 19 .
Figure 18.Dispersion and attenuation curves of Lamb waves from elastic plates coated with viscoelastic layer in water.(a) Dispersion curve; (b) attenuation curve.As shown in Figure19, the data of acoustic scattering echoes is measured from the side of steel; the time-angle figure is obtained by experiment for different frequencies of 800 kHz and 1.1 MHz.From Figures 15, 16 and 19, it is possible to clearly observe the highlights of mirror reflection echoes, edge wave echoes, and elastic scattering echoes.Among them, the edge waves are reflected in the diagram as two bright lines, indicating that the scattering echoes produced by the edges can be received by the transducer at most angles.The two edge waves intersect at an angle position in the pseudo-color image, and their intersection is the mirror reflection echo.The highlights of elastic scattering echoes generated by Lamb waves can also be distinctly observed in the image, appearing on the bright lines produced by the edge waves.Sensors 2023, 23, x FOR PEER REVIEW 19 of 23

Figure 19 .
Figure 19.Degree-time pseudo picture of acoustic echoes from the steel plate coated with a viscoelastic layer.(a) Frequency 800 kHz; (b) frequency 1.1 MHz.

Table 2 .
Highlights from steel plate 2.Another experiment is also conducted on a steel plate, named steel plate 2, whose size is 151.10mm133.20mm 1.87mm

Table 3 .
Highlights of the steel plate coated with a viscoelastic layer.Tables 1-3 present the theoretical calculations and experimental measurements of the critical angles for the highlights of elastic scattering echoes for steel plate 1, steel plate 2, and the coated steel plate.It can be observed from the corresponding results in the tables

Table 3 .
Highlights of the steel plate coated with a viscoelastic layer.