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

Authors: Yanhui Zhang; Qiang Gao
Journal: Shock and Vibration
 (2022)

---

## Abstract

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.


---

## Body

## 1. Introduction

Magneto-electro-elastic (MEE) materials, which can achieve the conversion between electrical and magnetic energy due to the coupled effect between electric and magnetic fields, have permeated every aspect of the modern technology. This feature promotes the wide application of the MEE in many fields of science and engineering [1–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–10] contributed many efforts 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-infinite 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 infinite MEE plate. Using the ordinary differential equation and stiffness matrix method, Ezzin et al. [16] investigated the propagation behavior of SH waves in laminated MEE plates, and the effects 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. Their 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. The 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–3].These 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. The 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. These 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.The 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. The combination of the PIM and the Wittrick–Williams (W-W) algorithm [28–31] has used extensively to study wave propagation in layered media [32–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. Then, 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. The performance of the method is verified by several numerical examples.


## 2. Statement of the Problem and Basic Equations

Consider a multilayered MEE structure consisting ofl layers with different material properties and layer thicknesses, as illustrated in Figure 1. The origin of the Cartesian coordinate system O−xyz is set on the surface, where the z-axis is along the depth direction. The structure extends infinitely in the x and y directions. The kth layer is bounded by the interfaces zk−1 and zk. The thickness for the kth layer is hk=zk+1−zk (k = 1, 2,…, l), and the total thickness is H=∑k=1lhk. 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.Figure 1
Schematic diagram for a multilayered MEE structure.For an anisotropic MEE solid, the constitutive equation in linear approximation can be expressed by equation [38].(1)σ¯D¯B¯=c−eT−hTeεαhαTμS¯E¯H¯,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(2)∇⋅σ¯=ρ∂2u¯∂t2,∇⋅D¯=0,∇⋅B¯=0,where u¯ is the displacement vector, t denotes time, ρ is the density, and ∇ is the nabla operator. The generalized geometric equations are described by(3)S¯=12∇u¯+∇u¯T,E¯=−∇φ¯,H¯=−∇ψ¯,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(4)Ω=εααTμ,and elastic coefficient matrix c are both positive definite, as shown in Appendix A.The plane-wave solution is considered as(5)q¯x,y,z,t=q^κx,κy,z,ωexpiκxx+κyy−ωt,where i=−1; q¯=u¯,v¯,w¯,φ¯,ψ¯T and q^=u¯,v¯,w¯,φ¯,ψ¯T 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(6)K22q″+K21−K12q′^−K11−ρω2Λq^=0,where q′^ and q″^ denote the first and second derivatives of q^ with respect to z; and(7)K22=c55c45c35e35h35c45c44c34e34h34c35c34c33e33h33e35e34e33−ε33−α33h35h34h33−α33−μ33,Λ=1000001000001000000000000,(8)K21=K12H=ikxc15c56c55e15h15c14c46c45e14h14c13c36c35e13h13e31e36e35−ε13−α31h31h36h35−α13−μ13+ikyc56c25c45e25h25c46c24c44e24h24c36c23c34e23h23e36e32e34−ε23−α32h36h32h34−α23−μ23,(9)K11=kx2c11c16c15e11h11c16c66c56e16h16c15c56c55e15h15e11e16e15−ε11−α11h11h16h15−α11−μ11+ky2c66c26c46e26h26c26c22c24e22h22c46c24c44e24h24e26e22e24−ε22−α22h26h22h24−α22−μ22+kxky2c16c12+c66c14+c56e21+e16h21+h16c12+c662c26c46+c25e26+e12h26+h12c14+c56c46+c252c45e25+e14h25+h14e21+e16e26+e12e25+e14−2ε12−α12+α21h21+h16h26+h12h25+h14−α12+α21−2μ12.Here, the superscript H denotes a Hermitian transpose. It can be seen thatK22 and K11 are both real symmetric matrices, and K12 and K21 are both pure imaginary matrices.


## 3. State-Space Formulation and Symplectic Transformation

### 3.1. State-Space Formulation

Introducing a dual vectorp^=τxz,τyz,σz,Dz,BzT satisfied by p^=K22q′^+K21q^ and defining the state vector v^=q^T,p^TT, from equation (6) the following equation of state can be given as(10)v′^=H^v^,H^=A^D^B^−A^H,where(11)A^=−K22−1K21,B^=K11−K12K22−1K21−ρω2Λ,D^=K22−1.It is easy to verify thatH^ is a Hamiltonian matrix, satisfied by H^H=JH^J, where J is a unit symplectic matrix defined by(12)J=0I5−I50.and In 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]. Then, combined to the W-W algorithm, all eigenfrequencies of waves in layered elastic media can be calculated [33, 34]. The 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. This 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 matrixc, from the form of K22 in equation (7) and the relation of D^=K22−1 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 K22 defined in equation (7) is not positive definite. Then, according to the relation of D^=K22−1, it is evident that D^ is not positive definite. The 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.


### 3.2. Symplectic Transformation

Introducing the following mathematical transformation,(13)v^=Tv,in which T and v can be expressed as(14)T=I303,203,303,202,302,202,3−I203,303,2I303,202,3I202,302,2,v=qp,q=u,v,w,Dz,BzT,p=τxz,τyz,σz,−φ,−ψT.It can be easily confirmed thatT is satisfied by TTJT=J, so it is a symplectic matrix [40]. Substituting equation (13) into equation (10) gives(15)v′=Hv,H=T−1H^T.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:(16)H=ADB−AH.After performing the aforementioned symplectic transformation, the new matrixD will be positive definite. The proof of positive definitiveness of matrix D is given in Appendix B.


## 3.1. State-Space Formulation

Introducing a dual vectorp^=τxz,τyz,σz,Dz,BzT satisfied by p^=K22q′^+K21q^ and defining the state vector v^=q^T,p^TT, from equation (6) the following equation of state can be given as(10)v′^=H^v^,H^=A^D^B^−A^H,where(11)A^=−K22−1K21,B^=K11−K12K22−1K21−ρω2Λ,D^=K22−1.It is easy to verify thatH^ is a Hamiltonian matrix, satisfied by H^H=JH^J, where J is a unit symplectic matrix defined by(12)J=0I5−I50.and In 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]. Then, combined to the W-W algorithm, all eigenfrequencies of waves in layered elastic media can be calculated [33, 34]. The 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. This 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 matrixc, from the form of K22 in equation (7) and the relation of D^=K22−1 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 K22 defined in equation (7) is not positive definite. Then, according to the relation of D^=K22−1, it is evident that D^ is not positive definite. The 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.


## 3.2. Symplectic Transformation

Introducing the following mathematical transformation,(13)v^=Tv,in which T and v can be expressed as(14)T=I303,203,303,202,302,202,3−I203,303,2I303,202,3I202,302,2,v=qp,q=u,v,w,Dz,BzT,p=τxz,τyz,σz,−φ,−ψT.It can be easily confirmed thatT is satisfied by TTJT=J, so it is a symplectic matrix [40]. Substituting equation (13) into equation (10) gives(15)v′=Hv,H=T−1H^T.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:(16)H=ADB−AH.After performing the aforementioned symplectic transformation, the new matrixD will be positive definite. The proof of positive definitiveness of matrix D is given in Appendix B.


## 4. 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. Those formulas may be considered as an extension of the purely elastic structure, and they can be found in previously published literature [33, 34, 36].The PIM [27] is first used to divide the kth layer into 2Nk sublayers with equal thickness γk=hk/2Nk, from which two arbitrary adjacent sublayers za,zb and zb,zc are selected. Then, the solution of equation (10) based on the mixed energy matrix form gives the following equations:(17)qbpa=Fa,bGa,b−Qa,bFa,bHqapb,(18)qcpb=Fb,cGb,c−Qb,cFb,cHqbpc,in which Fa,b, Ga,b, Qa,b and Fb,c, Gb,c, Qb,c are mixed energy matrices for the sublayer za,zb and zb,zc; qa, qb, qc and pa, pb, pc are the displacement and dual vectors at the interfaces za, zb and zc. Using equations (17) and (18) and eliminating qb and pb, the following equation for the sublayer za,zc is obtained as follows:(19)qcpa=Fa,cGa,c−Qa,cFa,cHqapc,where(20)Fa,c=Fb,cI+Ga,bQb,c−1Fa,b,Ga,c=Gb,c+Fb,cGa,b−1+Qb,cFb,cH,Qa,c=Qa,b+Fa,bHQb,c−1+Ga,b−1Fa,b.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. Then, a symbol composed of three letters is used to describe the boundary condition, for example, COO denotes the boundary condition of mechanically clamped, electrically open, and magnetically open. Assume thatJa,bCOO−COO and Jb,cCOO−COO are the eigenvalue counts for the sublayers za,zb and zb,zc with both ends mechanically clamped, electrically open, and magnetically open (COO-COO). Then, according to the W-W algorithm [28], the eigenvalue count for sublayer za,zc with both ends mechanically clamped, electrically open and magnetically open can be calculated by [29](21)Ja,cCOO−COO=Ja,bCOO−COO+Jb,cCOO−COO+sGa,b−1+Qb,c+Fb,cHGb,c−1Fb,c,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) Nk times, the eigenvalue count Jk−1,kCOO−COO and mixed energy matrices Fk−1,k, Gk−1,k and Qk−1,k for the kth layer can be obtained. Then, again repeated application of equations (20) and (21) l − 1 times, the eigenvalue count J0,lCOO−COO and global mixed energy matrices F0,l, G0,l and Q0,l for all l layers can be obtained.Note that the mixed energy matrices for a sublayer with thicknessγ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 large Nk can guarantee JγkCOO−COO=0 for a sublayer with thickness γk.


## 5. 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 Figure2) 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.Figure 2
Schematic diagram of the structure for five types of typical boundary conditions, where (a) Case (I) COO-COO, (b) Case II: FOO-FOO, (c) Case III: CSS-FOO, (d) Case IV: CSS-FSS, and (e) Case V: periodical structure.
(a)(b)(c)(d)(e)On the boundary surfacesz=z0 and z=zl, commonly encountered mechanical, electrical, and magnetical boundary conditions are considered as follows:Case I: The boundaries at both surfaces are mechanically clamped, electrically open, and magnetically open (COO-COO), that is,(22)uz0=vz0=wz0=Dzz0=Bzz0=0,uzl=vzl=wzl=Dzzl=Bzzl=0.The eigenvalue count for this case has been given in Section4, that is, JCOO−COO.Case II: The boundaries at both surfaces are mechanically free, electrically open, and magnetically open (FOO-FOO), that is,(23)τzxz0=τzyz0=σzz0=Dzz0=Bzz0=0,τzxzl=τzyzl=σzzl=Dzzl=Bzzl=0.According to the W-W algorithm, the eigenvalue count for this case is(24)J0,lFOO−FOO=J0,lCOO−COO+sG0,lMM−G0,lMEG0,lEE−1G0,lEM−1+sQ0,lMM+F0,lEMHG0,lEE−1F0,lEM.Case III: The boundary at the surfacez=z0 is mechanically clamped, electrically short and magnetically short, and the boundary at the surface z=zl is mechanically free, electrically open, and magnetically open (CSS-FOO), that is,(25)uz0=vz0=wz0=φz0=ψz0=0,τzxzl=τzyzl=σzzl=Dzzl=Bzzl=0.According to the W-W algorithm, the eigenvalue count for this case is(26)J0,lCSS−FOO=J0,lCOO−COO+sG0,lMM−G0,lMEG0,lEE−1G0,lEM−1+sQ0,lEE+F0,lEEHG0,lEE−1F0,lEE.Case IV: The boundary at the surfacez=z0 is mechanically clamped, electrically short and magnetically short, and the boundary at the surface z=zl is mechanically free, electrically short, and magnetically short (CSS-FSS), that is,(27)uz0=vz0=wz0=φz0=ψz0=0,τzxzl=τzyzl=σzzl=φzl=ψzl=0.According to the W-W algorithm, the eigenvalue count for this case is(28)J0,lCSS−FSS=J0,lCOO−COO+sG0,l−1+sQ0,lEE.Case V: For a periodically infinite multilayered structure, the state vectorsv0 and vl at both ends of a unit cell satisfy the Bloch’s theorem [42], that is,(29)vl=expiξHv0,whereξ is the Bloch wavenumber in periodically infinite multilayered structures. According to the W-W algorithm, the eigenvalue count for this case is(30)J∞,∞period=J0,lCOO−COO+sG0,l−1I5−e−iξHF0,l+I5−eiξHF0,lH−1Q0,l.So far, the eigenvalue count for the whole structure with a certain boundary condition has been obtained for a given wavenumberκ and trial frequency ω#. Then, 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.


## 6. Numerical Results and Discussion

Three numerical examples are considered in this section. The focus is on examining the numerical performance of the method. In Section6.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. The parameters of materials BaTiO3, CoFe2O4, Bimaterial A, and Bimaterial B used in examples are listed in Table 1 [8,13,43].Table 1
Material properties used in the examples:c (109 GPa), e (C/m2), h (N/Am), ε (10−9 C2/Nm2), μ (10−6 Ns2/C2), α (Ns2/C2), and ρ (kg/m3).
BaTiO3CoFe2O4BimaterialABc11166286217.6486.74c1277173125.3−8.25c1378170.5127.3127.15c14001.34−3.66c1500−2.870c1600−1.430c22166286222.02129.77c2378170.5125.44−7.42c2400−2.725.7c25000.490c2600−2.350c33162269.5216.25102.83c3400−1.139.92c3500−1.950c36001.620c444345.347.3138.81c45000.540c4600−1.680c554345.347.9468.81c56000.092.53c6644.556.547.829.01e11005.730.171e1200−1.35−0.152e1300−1.38−0.019e1400−0.020.067e1511.604.070e16002.040e2100−0.790e22003.320e2300−0.80e2411.604.090e2500−0.020.108e26003.54−0.095e31−4.40−1.580e32−4.40−1.560e3318.606.610e34002.030e35003.52−0.076e3600−0.020.067h1100401.950h1200140.190h130027.620h1400−75.050h15055064.47550h160032.240h210037.610h2200275.390h230015.950h240550151.13550h2500−750h2600130.880h310580.375.22580.3h320580.3161.87580.3h330699.7420.8699.7h340010.570h350018.310h3600−75.050ε1111.20.0860.392ε12ε12ε120000ε13000.350ε2211.20.085.640.398ε230000.86ε3312.60.09360.404μ115590190.255μ120000μ1300−106.80μ2255902975μ230000μ3310157190.2510α11000.050α120000α130000α22000.050α230000α33000.030ρ5800530055504900


### 6.1. A Multilayered Transversely Isotropic MEE Structure

A three-layered Sandwich plate with BaTiO3/CoFe2O4/BaTiO3 (B/F/B) configuration is considered in this example, and all three layers have equal thickness. The 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. The wave propagation direction is assumed to be along the x-axis, that is, θ=0o. To facilitate the comparison, the dimensionless wavenumber and dimensionless phase velocity are defined as κH and c¯p=ω/κcmax/ρmax; and cmax and ρmax are the maximum elastic constant and maximum density in the system.The dispersion curves for the B/F/B plate for four different types of boundary conditions are shown in Figure3, where subfigures (a), (b), (c), and (d) represent the results for the boundary conditions of FOO-FOO, COO-COO, CSS-FOO, and CSS-FSS. Generally, in the case of anisotropic MEE material [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 BaTiO3 and CoFe2O4 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) and 3(b), it can be seen that the dispersion curves have significant changes when the boundary condition changes from FOO-FOO to COO-COO. However, compared Figures 3(c) and 3(d), it can be observed that the dispersion curves are only slight variations when the boundary condition changes from CSS-FOO to CSS-FSS. This indicates that the effects of the mechanical boundary on dispersion curves are more pronounced than that of electrical and magnetical boundaries.Figure 3
Dispersion curves of guided waves in the B/F/B plate with the boundaries of (a) FOO-FOO, (b) COO-COO, (c) CSS-FOO, and (d) CSS-FSS.
(a)(b)(c)(d)Figure 4
Dispersion curves of guided waves in the three-layered plates with different stacking sequences for (a) the first mode, (b) the second mode, (c) the third mode, and (d) the fourth mode.
(a)(b)(c)(d)Next, the effects of varying stacking sequence to plate on dispersion curves are examined. The boundary conditions of the top and bottom surfaces are specified to FOO-FOO. The plates with stacking sequences B/B/B, B/F/B, F/B/F, and F/F/F are considered, respectively. The dispersion curves for the first four modes are shown in Figures4(a)–4(d), respectively, where the solid, dotted, dash-dotted, and dashed curves denotes the dispersion curves for stacking sequences B/B/B, B/F/B, F/B/F, and F/F/F, respectively. 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 along the z-axis direction, the semidiscretized equation can be obtained by(31)KΦ=ω2MΦ,in which stiffness matrix K and mass matrix M can be obtained by assembling the element counterparts. Then, equation (29) can be solved by using the eigs function in MATLAB for a given wavenumber κ. At κH=2, the dimensionless frequencies Ω=ωH/cmax/ρmax obtained from this method and the SAFE method discretized each layer with 0.001 m thickness sublayer are listed in Table 2. The results show very good agreement between the results obtained from the proposed method and the SAFE method.Table 2
AtκH=2, the first five dimensionless frequencies Ω of the three-layered plate with different stacking sequence for this method and SAFE method.
OrderB/B/BB/F/BF/B/FF/F/FPresentSAFEPresentSAFEPresentSAFEPresentSAFE10.72233000.72233000.54705030.54705030.59571870.59571870.56433710.564337221.03551391.03551390.83484240.83484240.88174770.88174770.88893750.888937431.74706411.74706411.41046871.41046871.45641481.45641471.44625001.446250041.90496021.90496021.47404821.47404821.57980771.57980781.53410311.534103152.54972402.54972411.92427581.92427592.12265112.12265121.98923711.9892372


### 6.2. A Multilayered Strongly Anisotropic MEE Structure

The 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. The structure consists of four layers with stacking sequences of materials BaTiO3, CoFe2O4, Bimaterial A, and Bimaterial B. Dispersion curves of the first 10 modes for θ=30o are shown in Figure 5, where (a), (b), (c), and (d) denote the results for the boundary conditions of FOO-FOO, COO-COO, CSS-FOO, and CSS-FSS, respectively. Comparing Figure 5(a) with Figure 5(b), we observe that the dispersion curves are significantly different when the boundary condition changes from FOO-FOO to COO-COO. However, comparing Figure 5(c) with Figure 5(d), we observe that the dispersion curves have only slight variations when the boundary condition changes from CSS-FOO to CSS-FSS. Therefore, 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.Figure 5
Dispersion curves of guided waves in the four-layered plate with the boundaries of (a) FOO-FOO, (b) COO-COO, (c) CSS-FOO, and (d) CSS-FSS.
(a)(b)(c)(d)For the dispersion diagram shown in Figure5(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.Figure 6
The mechanical displacements, electric and magnetic potentials of guided waves in the four-layered plate for the first three modes.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)Figure 7
The stress components, electrical and magnetical displacements of guided waves in the four-layered plate for the first three modes.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)Table 3
AtκH=2, the first five frequencies ω(rad·μs−1) of this method and the SAFE method for four types of boundary conditions.
OrderκHCOO-COOFOO-FOOCSS-FOOCSS-FSSPresentSAFEPresentSAFEPresentSAFEPresentSAFE10.12.19039172.19039180.00334430.00334491.14764031.14764041.14759821.147598212.23581912.23581920.08082960.08082961.20821221.20821221.20790611.207906152.34865112.34865110.29526570.29526571.37929811.37929811.37795871.3779587107.92618977.92618986.91305326.91305407.33915327.33915396.82912336.829124020.12.50261322.50261330.07550950.07550941.16257651.16257651.16252671.162526712.61855972.61855980.37752930.37752931.31612161.31612161.31432701.314327052.91697722.91697720.75493530.75493531.69731681.69731681.68934421.6893442108.09770348.09770377.32701137.32701167.72121337.72121347.70858427.708584330.14.41746894.41746920.12661320.12661322.52419622.52419632.51535212.515352114.36708994.36709010.63310440.63310442.48746382.48746392.47578532.475785454.27094684.27094701.25832861.25832862.42505172.42505182.40418212.4041822108.73323888.73323917.39726087.39726138.15180018.15180088.13164768.131648140.14.58668674.58668692.28576192.28576203.41021313.41021333.40918573.409185914.58134034.58134072.31591882.31591883.44126583.44126603.43964043.439640654.61598444.61598482.40508552.40508553.51467373.51467393.51273813.5127383109.46507569.46507727.82313157.82313168.48198038.48198058.47739968.477399950.14.85085754.85085802.40548182.40548183.60301223.60301243.60240393.602404115.00089205.00089242.48879222.48879233.71298263.71298283.70751363.707513955.27977645.27977682.72276252.72276263.98092283.98092303.97508323.97508341010.23627810.2362798.45652248.45652369.24658639.24658759.23690859.2369096


### 6.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. The unit cell consists of four layers with materials BaTiO3, CoFe2O4, Bimaterial A and Bimaterial B, and their thicknesses are the same. The incident angle θ is assumed to be 30o.The 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 Figure8, 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. The 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.Figure 8
Dispersion curves of guided waves in the periodic plates with 1 unit cell (dashed curves), 10 unit cells (solid curves), and 100 unit cells (dotted curves) under the four types of boundaries: (a) FOO-FOO, (b) COO-COO, (c) CSS-FOO, and (d) CSS-FSS.
(a)(b)(c)(d)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. Figure9 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. The 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.Figure 9
Dispersion curves (red solid lines) and bandgaps (shaded segments) for a periodically infinite multilayered structure, where (a)κH=0.1, (b) κH=1, (c) κH=5, and (d) κH=10.
(a)(b)(c)(d)


## 6.1. A Multilayered Transversely Isotropic MEE Structure

A three-layered Sandwich plate with BaTiO3/CoFe2O4/BaTiO3 (B/F/B) configuration is considered in this example, and all three layers have equal thickness. The 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. The wave propagation direction is assumed to be along the x-axis, that is, θ=0o. To facilitate the comparison, the dimensionless wavenumber and dimensionless phase velocity are defined as κH and c¯p=ω/κcmax/ρmax; and cmax and ρmax are the maximum elastic constant and maximum density in the system.The dispersion curves for the B/F/B plate for four different types of boundary conditions are shown in Figure3, where subfigures (a), (b), (c), and (d) represent the results for the boundary conditions of FOO-FOO, COO-COO, CSS-FOO, and CSS-FSS. Generally, in the case of anisotropic MEE material [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 BaTiO3 and CoFe2O4 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) and 3(b), it can be seen that the dispersion curves have significant changes when the boundary condition changes from FOO-FOO to COO-COO. However, compared Figures 3(c) and 3(d), it can be observed that the dispersion curves are only slight variations when the boundary condition changes from CSS-FOO to CSS-FSS. This indicates that the effects of the mechanical boundary on dispersion curves are more pronounced than that of electrical and magnetical boundaries.Figure 3
Dispersion curves of guided waves in the B/F/B plate with the boundaries of (a) FOO-FOO, (b) COO-COO, (c) CSS-FOO, and (d) CSS-FSS.
(a)(b)(c)(d)Figure 4
Dispersion curves of guided waves in the three-layered plates with different stacking sequences for (a) the first mode, (b) the second mode, (c) the third mode, and (d) the fourth mode.
(a)(b)(c)(d)Next, the effects of varying stacking sequence to plate on dispersion curves are examined. The boundary conditions of the top and bottom surfaces are specified to FOO-FOO. The plates with stacking sequences B/B/B, B/F/B, F/B/F, and F/F/F are considered, respectively. The dispersion curves for the first four modes are shown in Figures4(a)–4(d), respectively, where the solid, dotted, dash-dotted, and dashed curves denotes the dispersion curves for stacking sequences B/B/B, B/F/B, F/B/F, and F/F/F, respectively. 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 along the z-axis direction, the semidiscretized equation can be obtained by(31)KΦ=ω2MΦ,in which stiffness matrix K and mass matrix M can be obtained by assembling the element counterparts. Then, equation (29) can be solved by using the eigs function in MATLAB for a given wavenumber κ. At κH=2, the dimensionless frequencies Ω=ωH/cmax/ρmax obtained from this method and the SAFE method discretized each layer with 0.001 m thickness sublayer are listed in Table 2. The results show very good agreement between the results obtained from the proposed method and the SAFE method.Table 2
AtκH=2, the first five dimensionless frequencies Ω of the three-layered plate with different stacking sequence for this method and SAFE method.
OrderB/B/BB/F/BF/B/FF/F/FPresentSAFEPresentSAFEPresentSAFEPresentSAFE10.72233000.72233000.54705030.54705030.59571870.59571870.56433710.564337221.03551391.03551390.83484240.83484240.88174770.88174770.88893750.888937431.74706411.74706411.41046871.41046871.45641481.45641471.44625001.446250041.90496021.90496021.47404821.47404821.57980771.57980781.53410311.534103152.54972402.54972411.92427581.92427592.12265112.12265121.98923711.9892372


## 6.2. A Multilayered Strongly Anisotropic MEE Structure

The 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. The structure consists of four layers with stacking sequences of materials BaTiO3, CoFe2O4, Bimaterial A, and Bimaterial B. Dispersion curves of the first 10 modes for θ=30o are shown in Figure 5, where (a), (b), (c), and (d) denote the results for the boundary conditions of FOO-FOO, COO-COO, CSS-FOO, and CSS-FSS, respectively. Comparing Figure 5(a) with Figure 5(b), we observe that the dispersion curves are significantly different when the boundary condition changes from FOO-FOO to COO-COO. However, comparing Figure 5(c) with Figure 5(d), we observe that the dispersion curves have only slight variations when the boundary condition changes from CSS-FOO to CSS-FSS. Therefore, 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.Figure 5
Dispersion curves of guided waves in the four-layered plate with the boundaries of (a) FOO-FOO, (b) COO-COO, (c) CSS-FOO, and (d) CSS-FSS.
(a)(b)(c)(d)For the dispersion diagram shown in Figure5(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.Figure 6
The mechanical displacements, electric and magnetic potentials of guided waves in the four-layered plate for the first three modes.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)Figure 7
The stress components, electrical and magnetical displacements of guided waves in the four-layered plate for the first three modes.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)Table 3
AtκH=2, the first five frequencies ω(rad·μs−1) of this method and the SAFE method for four types of boundary conditions.
OrderκHCOO-COOFOO-FOOCSS-FOOCSS-FSSPresentSAFEPresentSAFEPresentSAFEPresentSAFE10.12.19039172.19039180.00334430.00334491.14764031.14764041.14759821.147598212.23581912.23581920.08082960.08082961.20821221.20821221.20790611.207906152.34865112.34865110.29526570.29526571.37929811.37929811.37795871.3779587107.92618977.92618986.91305326.91305407.33915327.33915396.82912336.829124020.12.50261322.50261330.07550950.07550941.16257651.16257651.16252671.162526712.61855972.61855980.37752930.37752931.31612161.31612161.31432701.314327052.91697722.91697720.75493530.75493531.69731681.69731681.68934421.6893442108.09770348.09770377.32701137.32701167.72121337.72121347.70858427.708584330.14.41746894.41746920.12661320.12661322.52419622.52419632.51535212.515352114.36708994.36709010.63310440.63310442.48746382.48746392.47578532.475785454.27094684.27094701.25832861.25832862.42505172.42505182.40418212.4041822108.73323888.73323917.39726087.39726138.15180018.15180088.13164768.131648140.14.58668674.58668692.28576192.28576203.41021313.41021333.40918573.409185914.58134034.58134072.31591882.31591883.44126583.44126603.43964043.439640654.61598444.61598482.40508552.40508553.51467373.51467393.51273813.5127383109.46507569.46507727.82313157.82313168.48198038.48198058.47739968.477399950.14.85085754.85085802.40548182.40548183.60301223.60301243.60240393.602404115.00089205.00089242.48879222.48879233.71298263.71298283.70751363.707513955.27977645.27977682.72276252.72276263.98092283.98092303.97508323.97508341010.23627810.2362798.45652248.45652369.24658639.24658759.23690859.2369096


## 6.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. The unit cell consists of four layers with materials BaTiO3, CoFe2O4, Bimaterial A and Bimaterial B, and their thicknesses are the same. The incident angle θ is assumed to be 30o.The 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 Figure8, 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. The 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.Figure 8
Dispersion curves of guided waves in the periodic plates with 1 unit cell (dashed curves), 10 unit cells (solid curves), and 100 unit cells (dotted curves) under the four types of boundaries: (a) FOO-FOO, (b) COO-COO, (c) CSS-FOO, and (d) CSS-FSS.
(a)(b)(c)(d)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. Figure9 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. The 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.Figure 9
Dispersion curves (red solid lines) and bandgaps (shaded segments) for a periodically infinite multilayered structure, where (a)κH=0.1, (b) κH=1, (c) κH=5, and (d) κH=10.
(a)(b)(c)(d)


## 7. 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. The feasibility of the W-W algorithm after performing the symplectic transformation is demonstrated by a strictly theoretical analysis. The 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.The 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. The 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.



---

*Source: 1346719-2022-06-06.xml*