A Generalized Method for Dispersion Analysis of Guided Waves in Multilayered Anisotropic Magneto-Electro-Elastic Structures

Based on the symplectic structure of the Hamiltonian matrix, the precise integration method (PIM), and the Wittrick–Williams (W-W) algorithm, a generalized method for computing the dispersion curves of guided waves in multilayered anisotropic magneto-electro-elastic (MEE) structures for different types of mechanical, electrical, and magnetical boundaries is developed. A strictly theoretical analysis shows that the W-W algorithm cannot be applied directly to the MEE structure. This is because a block of the Hamiltonian matrix is not positive definite for MEE structures so that the eigenvalue count of the sublayer is not zero when the divided sublayer is sufficiently thin. To overcome this difficulty, based on the symplectic structure of the Hamiltonian matrix, a symplectic transformation is introduced to ensure that the W-W algorithm can be applied conveniently to solve wave propagation problems in multilayered anisotropic MEE structures. The application of the PIM based on the mixed energy matrix to solve the wave equation can ensure the stability and efficiency of the method, and all eigenfrequencies are found without the possibility of any being missed using the W-W algorithm. This research provides the necessary insight to apply the W-W algorithm in wave propagation and vibration problems of MEE structures.


Introduction
Magneto-electro-elastic (MEE) materials, which can achieve the conversion between electrical and magnetic energy due to the coupled e ect between electric and magnetic elds, have permeated every aspect of the modern technology. is feature promotes the wide application of the MEE in many elds of science and engineering [1][2][3], for example, sensors [4], smart devices [5], and nondestructive evaluation [6], which are closely related to the knowledge of wave propagation. Various techniques related to the wave propagation in MEE structures were developed by numerous scholars. Pan et al. [7][8][9][10] contributed many e orts to the wave propagation in MEE structures, and they derived a series of analytical expressions for various MEE structures. Ezzin et al. [11] investigated Love waves in a transversely isotropic piezoelectric layer bonded to an MEE semi-in nite space and derived explicit dispersion equations. Li et al. [12] investigated propagation behavior of the Bleusteine-Gulyaev waves in a functionally graded transversely isotropic MEE half-space and derived analytically the dispersion equations under electromagnetically open and shorted conditions. An analytical treatment based on the transfer matrix was presented by Chen et al. [13] for the propagation of harmonic waves in multilayered MEE plates. An exact solution for shear horizontal (SH) waves propagating in a transversely isotropic MEE was derived by Nie et al. [14].
Apart from analytical solutions, abundant numerical methods were also developed for wave propagation in layered MEE structures. Wu et al. [15] explored the dispersive behavior of the symmetric and antisymmetric Lamb waves in an in nite MEE plate. Using the ordinary di erential equation and sti ness matrix method, Ezzin et al. [16] investigated the propagation behavior of SH waves in laminated MEE plates, and the e ects of thickness ratio on phase velocity and group velocity were discussed. Based on Hamilton's principle, Xiao et al. [17,18] applied Chebyshev spectral element method to investigate the dispersion characteristics of guided waves in a functionally graded MEE plate and in a multilayered MEE curved panel, respectively. Chen et al. [19] employed the reverberation matrix method to compute the dispersion curves of functionally graded materials. eir method adopted a homogenous assumption of material parameters in each thin layer because the graded layer was divided into a sufficiently large number of thin layers. Yu et al. [20] employed Legendre orthogonal polynomial series method to study the dispersion behavior of guided waves in a layered MEE plate. Matar et al. [21] used the combined Legendre-Laguerre polynomial expansion method to compute dispersion curves and mode shapes in layered MEE composites. In addition, the study of periodic MEE structures has attracted wide attention in the past decades. Pang et al. [22,23] applied the transfer matrix method and the stiffness matrix method to investigate SH bulk/surface waves propagating in periodically layered infinite/semi-infinite MEE composites. Wang et al. [24] studied the wave propagation in periodic composites consisting of MEE plates using the plane-wave expansion method. Liu et al. [25] studied dispersion behavior and transmission coefficients of SH wave in a periodically layered piezoelectric structure and discussed the feature of bandgaps in periodic structures. Chen et al. [26] investigated dispersions and band structures of elastic waves in nanoscale periodic MEE structures based on the nonlocal theory. e localization factors and dispersion curves were computed using the transfer matrix method. Recently, the band structure and evanescent behavior of Bloch waves in periodic MEE structures were investigated by applying extended plane-wave expansion method and Bloch's theory [1][2][3].
ese numerical methods can be broadly classified into two groups. One group is based on analytical models, such as the transfer matrix method. In principle, these methods are exact because mathematical expressions are free from approximations, but in numerical calculation here exist many difficulties that are associated with the numerical overflow and the loss of precision at high frequencies. Meanwhile, a disadvantage is that these methods involve a challenging problem of seeking the roots from transcendental eigenequations, which means that its solution is usually obtained by the intensive search technique. e other group is based on discrete models, such as the finite-element method. For these methods, the substructuring technique is usually used to divide each layer into a large number of thin layers, then the displacements are approximate via interpolation polynomials. ese discrete methods avoid seeking the roots from the transcendental eigenequation by solving a quadratic eigenvalue problem. In spite of this advantage, it should be highlighted that, due to the need for the discretization, these methods give rise to more degrees of freedom of system than previous methods. e precise integration method (PIM) was proposed by Zhong [27] to solve accurately the sets of first-order ordinary differential equations with specified two-point boundary value conditions for space domain problems, or with specified initial value conditions for time-domain problems. e combination of the PIM and the Wittrick-Williams (W-W) algorithm [28][29][30][31] has used extensively to study wave propagation in layered media [32][33][34][35][36]. Numerous results shown that the method for solving wave propagation problem not only can avoid seeking the roots from the nonlinear transcendental eigenequation directly but also can guarantee that the calculations have good stability and high accuracy. Later, numerous researches show that the method based on the mixed energy matrix was more stable compared to those based on the transfer and stiffness matrices [36]. A key process for the application of the W-W algorithm is to determine the eigenvalue count J of the structure, and it is usually difficult to directly solve J. To this end, a substructure technique is usually used to divide each layer into sufficient number of sublayers, so that the eigenvalue count of the sublayer with both ends clamped is zero. en, by gradually condensing sublayers, the eigenvalue count of the whole structure is obtained. However, for a piezoelectric structure, the W-W algorithm is no longer suitable because the eigenvalue count of sublayer with both ends mechanically clamped, electrically short and magnetically short is zero even if each layer is divided into sufficient number of sublayers. For this reason, a symplectic transformation was adopted to deal with guided wave propagation in multilayered anisotropic piezoelectric structures [37]. After performing symplectic transformation, the W-W algorithm combined to the PIM could be applied to compute the eigenfrequencies of waves in the piezoelectric structure.
Obviously, similar problems exist in MEE structures as well when the W-W algorithm is applied to compute the wave dispersion. Compared to piezoelectric cases, solutions to the wave propagation problem in MEE structure are much more complicated because of its complex boundary conditions and high-dimensional governing equations. In this paper, based on symplectic transformation, a generalized method for dispersion analysis of waves in multilayered anisotropic MEE structures is developed. e performance of the method is verified by several numerical examples.

Statement of the Problem and Basic Equations
Consider a multilayered MEE structure consisting of l layers with different material properties and layer thicknesses, as illustrated in Figure 1. e origin of the Cartesian coordinate system O − xyz is set on the surface, where the z-axis is along the depth direction. e structure extends infinitely in the x and y directions. e kth layer is bounded by the interfaces z k− 1 and z k . e thickness for the kth layer is h k � z k+1 − z k (k � 1, 2,. . ., l), and the total thickness is H � l k�1 h k . Without the loss of generality, it is assumed to the wave propagating along the x-y plane and the incident angle θ measured from the positive x-axis in the clockwise direction.
For an anisotropic MEE solid, the constitutive equation in linear approximation can be expressed by equation [38].
where σ and S are the stress and strain tensors; D and B are electric displacement and magnetic induction tensors; E and H are the tensors of the electric and magnetic fields; c, e, and h are the elastic, piezoelectric, and piezomagnetic coefficient tensors; ε, μ, and α dielectric permittivity, magnetic permeability, and magneto-electric coefficient tensors; and the superscript T represents the transpose. In the absence of the body force, electric charge, and magnetic charge, the mechanical, electric, and magnetic governing equations can be expressed by where u is the displacement vector, t denotes time, ρ is the density, and ∇ is the nabla operator. e generalized geometric equations are described by where φ and ψ are the electric and magnetic potentials. Note that the detailed forms of the physical quantities and material constants for elastic, electric, and magnetic fields are presented in Appendix A. To ensure stability of the system, the total internal energy is more than zero [39], such that the magneto-electric matrix and elastic coefficient matrix c are both positive definite, as shown in Appendix A. e plane-wave solution is considered as are the generalized displacement vectors in the time-space and frequency-wavenumber fields, respectively; ω represents the wave circular frequency; κ x and κ y are the two components of the wave vector in the x and y directions, described separately as κ x � κ cos θ and κ y � κ sin θ; and κ is the magnitude of the wave vector along the propagation direction. Substituting equations (1) and (3) into equation (2) and then using equation (5), we obtain where q ′ and q ″ denote the first and second derivatives of q with respect to z; and

Shock and Vibration
Here, the superscript H denotes a Hermitian transpose. It can be seen that K 22 and K 11 are both real symmetric matrices, and K 12 and K 21 are both pure imaginary matrices. (6) the following equation of state can be given as

State-Space Formulation and Symplectic Transformation
where It is easy to verify that H is a Hamiltonian matrix, satisfied by H H � JHJ, where J is a unit symplectic matrix defined by and I n denotes a n × n identify matrix.
In the aforementioned paragraph, the governing equation of the wave motion was converted into a first-order ordinary differential equations, which can be accurately solved using the PIM [27]. en, combined to the W-W algorithm, all eigenfrequencies of waves in layered elastic media can be calculated [33,34]. e PIM is used to ensure the precision and stability of computation, and the W-W algorithm is used to ensure that all eigenfrequencies are found without the possibility of any being missed. However, the method could not be applied directly to wave propagation problems in MEE structures. is is due primarily to the fact that a block matrix D of the Hamiltonian matrix in equation (10) is not positive definite.
For purely elastic structures, according to the positive definiteness of the elastic coefficient matrix c, from the form of K 22 in equation (7) and the relation of D � K − 1 22 in equation (11), we conclude that D is positive definite. However, for MEE structures, matrices c and Ω are both positive definite (see Appendix A), so K 22 defined in equation (7) is not positive definite. en, according to the relation of D � K − 1 22 , it is evident that D is not positive definite. e computational formulas of the W-W algorithm require that D must be a positive definite matrix. In the following section, a symplectic transformation is introduced to overcome this difficulty.

Symplectic Transformation.
Introducing the following mathematical transformation, v � Tv, (13) in which T and v can be expressed as It can be easily confirmed that T is satisfied by T T JT � J, so it is a symplectic matrix [40]. Substituting equation (13) into equation (10) gives According to the properties of the Hamiltonian matrix and symplectic theory [40], if T is a symplectic matrix, H is also a Hamiltonian matrix, which can be written in the following form: After performing the aforementioned symplectic transformation, the new matrix D will be positive definite.
e proof of positive definitiveness of matrix D is given in Appendix B.

PIM and W-W Algorithm for Multilayered MEE Structures
After performing the symplectic transformation, this section will be devoted to the derivation of formulas for calculating wave propagation problems in multilayered MEE structures based on the PIM and the W-W algorithm. ose formulas may be considered as an extension of the purely elastic structure, and they can be found in previously published literature [33,34,36]. e PIM [27] is first used to divide the kth layer into 2 N k sublayers with equal thickness c k � h k /2 N k , from which two arbitrary adjacent sublayers [z a , z b ] and [z b , z c ] are selected.
en, the solution of equation (10) based on the mixed energy matrix form gives the following equations:  (17) and (18) and eliminating q b and p b , the following equation for the sublayer [z a , z c ] is obtained as follows:

Shock and Vibration
where For the sake of brevity, assuming that F and C represent the boundary conditions of mechanically free and clamped, and that O and S represent the boundary conditions of electrically or magnetically open and short. en, a symbol composed of three letters is used to describe the boundary condition, for example, COO denotes the boundary con- en, according to the W-W algorithm [28], the eigenvalue count for sublayer [z a , z c ] with both ends mechanically clamped, electrically open and magnetically open can be calculated by [29] where s # { } is the sign count of the matrix #, which is defined as the number of the negative eigenvalues of the matrix # and is equal to the number of changes in sign of the Sturm sequence [41].
Repeated application of equations (20) and (21) N k times, the eigenvalue count J COO− COO k− 1,k and mixed energy matrices F k− 1,k , G k− 1,k and Q k− 1,k for the kth layer can be obtained. en, again repeated application of equations (20) and (21) l − 1 times, the eigenvalue count J COO− COO 0,l and global mixed energy matrices F 0,l , G 0,l and Q 0,l for all l layers can be obtained.
Note that the mixed energy matrices for a sublayer with thickness c k can be calculated using the Taylor series, and the computational formula can be found from literature [33,36]. Moreover, since the matrix D is positive definite, sufficiently

Computing Eigenfrequencies for Different Boundary Conditions
In the aforementioned section, the procedure for computing the eigenvalue count for a multilayered structure with both ends mechanically clamped, electrically open, and magnetically open has been given. In this section, five types of mechanical, electrical, and magnetical boundary conditions are considered (see Figure 2) and under these boundary conditions, the computational formulas of the eigenvalue count of the multilayered structure are also given. Finally, all eigenfrequencies can be calculated based on the concept of the eigenvalue count.
On the boundary surfaces z � z 0 and z � z l , commonly encountered mechanical, electrical, and magnetical boundary conditions are considered as follows: Case I: e boundaries at both surfaces are mechanically clamped, electrically open, and magnetically open (COO-COO), that is, e eigenvalue count for this case has been given in Section 4, that is, J COO− COO . Case II: e boundaries at both surfaces are mechanically free, electrically open, and magnetically open (FOO-FOO), that is, According to the W-W algorithm, the eigenvalue count for this case is Shock and Vibration 5 Case III: e boundary at the surface z � z 0 is mechanically clamped, electrically short and magnetically short, and the boundary at the surface z � z l is mechanically free, electrically open, and magnetically open (CSS-FOO), that is, According to the W-W algorithm, the eigenvalue count for this case is Case IV: e boundary at the surface z � z 0 is mechanically clamped, electrically short and magnetically short, and the boundary at the surface z � z l is mechanically free, electrically short, and magnetically short (CSS-FSS), that is, According to the W-W algorithm, the eigenvalue count for this case is Case V: For a periodically infinite multilayered structure, the state vectors v 0 and v l at both ends of a unit cell satisfy the Bloch's theorem [42], that is, where ξ is the Bloch wavenumber in periodically infinite multilayered structures. According to the W-W algorithm, the eigenvalue count for this case is So far, the eigenvalue count for the whole structure with a certain boundary condition has been obtained for a given wavenumber κ and trial frequency ω # . en, based on the concept of the eigenvalue count, all eigenfrequencies can be computed using the bisection method [28,36]. Note that the general anisotropic case are derived in this paper, but the other special cases (e.g., antiplane and in-plane wave problems) can be easily recovered by removing some items of wave equation and reducing the size of matrices.

Numerical Results and Discussion
ree numerical examples are considered in this section. e focus is on examining the numerical performance of the method. In Section 6.1, an example for a three-layered MEE plate is first presented to verify the correctness of this method, which involves the comparison with previously published results. In Section 6.2, this method is used to compute the dispersion curves for a strongly anisotropic four-layered MEE plate. In Section 6.3, an example is presented to compute the dispersion curves of a periodically multilayered finite MEE structure and bandgaps of a periodically infinite multilayered MEE structure. e parameters of materials BaTiO 3 , CoFe 2 O 4 , Bimaterial A, and Bimaterial B used in examples are listed in Table 1 [8,13,43].

A Multilayered Transversely Isotropic MEE Structure.
A three-layered Sandwich plate with BaTiO 3 /CoFe 2 O 4 / BaTiO 3 (B/F/B) configuration is considered in this example, and all three layers have equal thickness. e objective of this example is to verify the correctness of this method compared to the published results, and then to examine the effects of varying boundary condition and stacking sequence to the plate on dispersion curves. e wave propagation direction is assumed to be along the x-axis, that is, θ � 0 o . To facilitate  [44,45], all waves with respect to the mechanical, electrical and magnetical components are coupled and in order to identify them one needs to solve a tenth order polynomial characteristic equation. For special material symmetry directions, such as isotropic or transversely isotropic materials, the wave motion can be decoupled into two sets of independent equations. One of the equations is a second-order polynomial equation corresponding to the antiplane wave with respect to the component v. Another equation is an eighth-order polynomial equation corresponding to the in-plane wave with respect to the component u, w, φ, and ψ. Since BaTiO 3 and CoFe 2 O 4 are both transversely isotropic materials, the wave motion can be decoupled into two sets of independent equations, corresponding to antiplane wave with respect to the component v and in-plane wave with respect to the components u, w, φ, and ψ. Figure 3(a) is in good agreement with the results given in Figures 4(c) of Ref. [13], which confirms the correctness of this method. From Figures 3(a) Figure 4(a) shows that for the first mode, the dispersion curves for the four different stacking sequences are close to each other. However, from Figures 4(b)-4(d), it is clear that for the second, third, and fourth modes, the dispersion curves are quite different for the four different stacking sequences.
To further demonstrate the performance of the proposed methods, the results obtained from the proposed methods are compared with those obtained from the semianalytical finite-element (SAFE) method [46]. Using the linear element   Shock and Vibration along the z-axis direction, the semidiscretized equation can be obtained by in which stiffness matrix K and mass matrix M can be obtained by assembling the element counterparts. en, equation (29) can be solved by using the eigs function in MATLAB for a given wavenumber κ. At κH � 2, the dimensionless frequencies Ω � ωH/ �������� c max /ρ max obtained from this method and the SAFE method discretized each layer with 0.001 m thickness sublayer are listed in Table 2. e results show very good agreement between the results obtained from the proposed method and the SAFE method.

A Multilayered Strongly Anisotropic MEE Structure.
e aforementioned example focuses on the transversely isotropic material. To confirm the feasibility of this method for the application of more generally anisotropic problem, a strongly anisotropic multilayered MEE structure is further considered in this section. e structure consists of four layers with stacking sequences of materials BaTiO 3  erefore, this example obtains similar conclusion to Section 6.1, that is, the effects of mechanical boundary are more significant than that of electrical and magnetical on dispersion curves.

Shock and Vibration
For the dispersion diagram shown inFigure 5(a), at κH � 1, the first three eigenfrequencies are ω 1 � 0.295266 rad·μs − 1 , ω 2 � 0.754935 rad·μs − 1 and ω 2 � 1.258329 rad·μs − 1 , and their corresponding eigenfunctions are shown in Figures 6 and 7. It can be seen from Figures 6 and 7 that all of mechanical, electrical, and magnetical components are complex functions, and their real and imaginary parts are separately plotted in (a)-(e) and in (f)-(g). From Figures 6 and 7, we observe that all these components u, v, w, φ, and ψ are not identically zero, which indicates that the wave motion is coupled with regard to all components. Furthermore, at κH � 0.1, 1, 5, and 10, the eigenfrequencies in the first five modes for four types of boundary conditions are listed in Table 3, which are provided for examination and reference by other researchers. In addition, to confirm the accuracy of this method, at κH � 0.1, 1, 5, and 10, the first five frequencies of this method and the SAFE method discretized each layer with 0.001 m thickness sublayer for four types of boundary conditions are listed in Table 3.

A Periodically Multilayered Anisotropic MEE Structure.
In this example, the results for both periodically multilayered finite and infinite MEE structures obtained from the presented method are presented. e unit cell consists of four layers with materials BaTiO 3 , CoFe 2 O 4 , Bimaterial A and Bimaterial B, and their thicknesses are the same. e incident angle θ is assumed to be 30 o . e dispersion curves of the first 10 modes for the periodically finite MEE structure for different types of boundary conditions and different number of unit cells are shown in Figure 8, where the dashed, solid, and dotted curves denote the phase velocities for the periodic multilayered plate composed of 1, 10, and 100 unit cells, respectively. e subfigures (a)-(d) represent the dispersion curves for the boundary conditions of FOO-FOO, COO-COO, CSS-FOO, and CSS-FSS. From Figure 8, we observe that the phase velocities for the structure consisting of 1 unit cell are significantly different from that consisting of 10 unit cells. However, the phase velocities for the structure consisting of 10 and 100 unit cells have small changes.
Finally, we investigate a periodic multilayered structure arranged by an infinite sequence of unit cell, and the unit cell is identified to that given in the finite periodic structure. Wave propagation in a periodical structure can exhibit a characteristic feature of band structure, known as passbands or Bragg bandgaps, within which the waves can propagate to infinity. In complementary frequency bands, known as stopbands or locally resonant bandgaps, the waves are effectively attenuated. Figure 9 shows the dispersion curves and band structures of the periodically infinite multilayered    structure in the first 10 modes, where (a), (b), (c), and (d) denote the results for κH � 0.1, 1, 5 and 10, respectively. e bandgaps are closely related to κH. When κH ≠ 0, the first bandgap starts from zero frequency. Only when κH � 0, the first bandgap does not start from zero frequency. From Figure 9, it can be observed that as κH increases, the width of the first stopband becomes large and the number of stopbands exhibits an increasing trend.

Conclusions
By combining the PIM with the W-W algorithm, based on the symplectic structure of the Hamiltonian matrix, a generalized method for calculating dispersion curves of waves in multilayered anisotropic MEE structures was proposed. To facilitate the application of the W-W algorithm in MEE structures, based on the symplectic structure of the Hamiltonian matrix, a symplectic transformation is introduced. e feasibility of the W-W algorithm after performing the symplectic transformation is demonstrated by a strictly theoretical analysis. e method combines the advantages of the PIM with those of the W-W algorithm so that the computation is accurate and stable, and no eigenfrequency is missed. Several typical examples are presented to verify the numerical performance of this method. In examples, the dispersion curves of waves in multilayered anisotropic MEE structures are calculated by using this method. It was shown that this method is highly accurate by comparing with the SAFE method. e proposed method can provide the basis of theory for the application of the PIM and the W-W algorithm to solve the wave propagation and dynamic problems of MEE structure. e method is not restricted to the specific material properties, the layer number of unit cells, and the complexity of the unit cell. Based on the idea of the proposed method, it can be easily applied to more complex problems. For example, by combining with the finite-element method, the proposed method can be extended to solve dispersions or bandgaps of finite, infinite, and semi-infinite anisotropic MEE phononic crystals.

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare no conflicts of interest.