The Computation of Complex Dispersion and Properties of Evanescent Lamb Wave in Functionally Graded Piezoelectric-Piezomagnetic Plates

Functionally graded piezoelectric-piezomagnetic material (FGPPM), with a gradual variation of the material properties in the desired direction(s), can improve the conversion of energy among mechanical, electric, and magnetic fields. Full dispersion relations and wave mode shapes are vital to understanding dynamic behaviors of structures made of FGPPM. In this paper, an analytic method based on polynomial expansions is proposed to investigate the complex-valued dispersion and the evanescent Lamb wave in FGPPM plates. Comparisons with other related studies are conducted to validate the correctness of the presented method. Characteristics of the guided wave, including propagating modes and evanescent modes, in various FGPPM plates are studied, and three-dimensional full dispersion and attenuation curves are plotted to gain a deeper insight into the nature of the evanescent wave. The influences of the gradient variation on the dispersion and the magneto-electromechanical coupling factor are illustrated. The displacement amplitude and electric potential and magnetic potential distributions are also discussed in detail. The obtained numerical results could be useful to design and optimize different sensors and transducers made of smart piezoelectric and piezomagnetic materials with high performance by adjusting the gradient property.


Introduction
Due to the excellent coupling behavior among mechanical, electric, and magnetic fields, piezoelectric-piezomagnetic composites (or magneto-electro-elastic material) composed of piezoelectric and piezomagnetic phases have been increasingly applied to different engineering structures, especially to the smart or intelligent systems as intelligent sensors, damage detectors, etc. [1]. It is found that the smart structures made of functionally graded materials (FGM) possess a better structural performance than traditional composite materials. The concept of FGM has been extended to the development of new piezoelectric-piezomagnetic materials appointed functionally graded piezoelectric-piezomagnetic materials, which can realize the smooth transition of the physical constitutive parameters of the piezoelectric and piezomagnetic materials. FGPPMs have been used in some devices to improve their efficiency and other features. Many applications are closely connected with the vibration and wave propagation of FGM and FGPPM [2][3][4][5]. Dispersion relations and wave mode shapes are very important for understanding dynamic behaviors of structures. Wave propagation features in FGPPM plates could be also useful in designing and optimizing the high-accuracy sensors and transducers [6,7].
With the remarkable achievements in fabrication of FGPPM during the decades, many investigators have turned attention to the study of wave propagation in such materials. Wang and Rokhlin [8] Materials 2018, 11, x FOR PEER REVIEW 3 of 16 direction with a thickness h, occupies the region 0 ≤ z ≤ h, as shown in Figure 1. The wave propagates along the x direction, and the upper and bottom surface of the plate are traction free. For the piezoelectric-piezomagnetic medium, the governing field equations (Equation (1)) and the generalized constitutive relations (Equation (2)) can be expressed as [12]: Generalized geometric equations under a rectangular coordinates system are in the above Equations (1)-(3), ij σ is the stress tensor, i D is the electric displacement and i B is the magnetic induction , , , , ijkl kij ijk ik ik C e q g ∈ and ik μ are the elastic, piezoelectric, piezomagnetic, dielectric, magnetic, and magnetoelectric parameters of the FGPPM, respectively, while all of them, including the mass density ρ, are functions of z. According to Einstein summation convention, where i, j, k and l = 1, 2, 3 corresponding to x, y, z directions, respectively.
where Ui (i = x, z) represents the amplitude of the displacements in the ith directions, X and Y represent the amplitude of electric potential and magnetic potential, respectively. k is the wave number, ω is the angular frequency, and i is the imaginary number.
Since the material properties change gradually with thickness and are the functions of z, they can be fitted into the following form: where f (f = ρ, C, e, ∈, q, g and μ) denotes material parameters, l is the order number, f (l) is the coefficient. For homogeneous material, , and when l > 0, f(l) is zero.
The following boundary and continuous conditions should be satisfied as follows. For the traction-free boundary condition, it requires that  For the piezoelectric-piezomagnetic medium, the governing field equations (Equation (1)) and the generalized constitutive relations (Equation (2)) can be expressed as [12]: Generalized geometric equations under a rectangular coordinates system are in the above Equations (1)-(3), σ ij is the stress tensor, D i is the electric displacement and B i is the magnetic induction C ijkl , e kij , q ijk , ∈ ik , g ik and µ ik are the elastic, piezoelectric, piezomagnetic, dielectric, magnetic, and magnetoelectric parameters of the FGPPM, respectively, while all of them, including the mass density ρ, are functions of z. According to Einstein summation convention, where i, j, k and l = 1, 2, 3 corresponding to x, y, z directions, respectively. ε kl , E k and H k are the strain tensor, the electric field, and magnetic field, respectively. u i (i = x, y, z) denotes the mechanical displacement component in the ith direction. ϕ and ψ are the electric potential and magnetic potential. Comma in subscripts and superposed dot denote spatial and time derivatives, respectively.
For Lamb waves propagating along the x direction, the displacement components, electric potential and magnetic potential can be expressed as where U i (i = x, z) represents the amplitude of the displacements in the ith directions, X and Y represent the amplitude of electric potential and magnetic potential, respectively. k is the wave number, ω is the angular frequency, and i is the imaginary number.
Since the material properties change gradually with thickness and are the functions of z, they can be fitted into the following form: where f (f = ρ, C, e, ∈, q, g and µ) denotes material parameters, l is the order number, f (l) is the coefficient.
Then take the traction-free and electrical and magnetic open-circuit boundary conditions as an example. Considering the boundary of material, the position-dependent material parameters are given by: (6) where π(z) is a rectangular window function defined by π(z) = 1, 0 ≤ z ≤ h 0, elsewhere , whose derivative is a Dirac's delta function, δ(z − h) − δ(z). Then the boundary conditions can be automatically incorporated in the constitutive relations [25].
To reduce the number of resolving equations, we substitute Equations (3)-(6) into Equation (2) with following substitution into Equation (1). Consequently, the governing differential equations in terms of the displacement, electric potential and magnetic potential components can be obtained. Here, the case of an orthotropic FGPPM plate with the z direction polarization is given: where U and W respectively represent the amplitude of vibration in the x and z directions. The superscript (') is the derivative with respect to z.
The four amplitudes can be expanded into Legendre orthogonal polynomial series as: where p α m (α = 1, 2, 3, 4) are the expansion coefficients, Q m (r) are an orthonormal set of polynomials in the interval [0,h].
where P m is the Legendre polynomial of order m. Substituting Equations (8) and (9) into Equation (7), then multiplying both sides of the modified Equation (7) by the complex conjugate Q * j (z) with j running from 0 to M, integrating over z from 0 to h, taking advantage of the orthonormality of the polynomial, yields: or is abbreviated as ∂z n dz. The objective is to find wave numbers k that satisfy the Equation (11). It is simple and useful for propagating wave, by specifying real k and then solving for ω. But if interest is the evanescent wave, the approach is useless because k is complex and the solving of Equation (11) involves a multivariable search. In order to overcome this difficulty, we develop a new solution procedure as shown below.
We introduce two new vectors: Substitution Equation (12) into Equation (11), and then multiplying both sides of the modified Equation (11) by inverse matrix A −1 , yields Combining Equation (13) and the above vector q = k · p, we obtain where I is the identity matrix and Z is a zero matrix.
If we define R = [pq] T , then Equation (14) can be written as Up to this stage, the problem is reduced to a typical eigenvalue problem, which can be easily solved using an eigensolver routine that yields the complex eigenvalues k. All the developments performed in this paper were implemented in Mathematica software (version 8.0, Wolfram company, Champaign, IL, USA). The calculation technique in the short-circuit case is similar to that which is used in the open-circuit case. The deduction process is not shown to save space.

Numerical Results and Discussion
Based on the previous formulations, the computer program in terms of the presented method has been written using Mathematica software to calculate the dispersion and phase velocity curves for the FGPPM plate composed of CoFe 2 O 4 (top) and Ba 2 TiO 3 (bottom), h = 1 mm. The material parameters are from literature [26] and are listed in Table 1. Table 1. Material parameters of two piezoelectric-piezomagnetic materials.

Materials
Property Units: C ij (10 9 N/m 2 ), ∈ ij (10 −10 F/m), e ij (C/m 2 ), q ij (N/Am), µ ij (10 −6 Ns 2 /C 2 ), ρ (10 3 kg/m 3 ). We use the Voigt-type model, as described in the literature [27], to calculate the effective material property of the FGPPM plate: where F B and F C respectively represent the material property of the Ba 2 TiO 3 and CoFe 2 O 4 materials, and V B and V C are volume fraction. Equation (16) can be rewritten as Similar to Equation (5), V C (z) can be expressed as a power expansion, Here we consider four different gradient fields, V C (z) = (z/h) n , n = 1, 2 and 3, namely linear, quadratic and cubic graded fields, and sinusoidal graded field V C (z) = sin(0.5πz/h).

Approach Validation and Convergence of the Problem
To check the validity and the efficiency of our approach, we make a comparison between our results and the literature results. Because there is no investigation on the evanescent waves in FGPPM so far, we compute the full spectrum of Lamb wave in a steel plate and make a comparison with the available results in literature [17] from a spectral collocation method. The calculating parameters are ρ = 7932 kg/m 3 , C 11 = 281.757 GPa, C 12 = 113.161 GPa, C 44 = 84.298 GPa, and h = 10 mm. The non-dimensional frequency and wave number are defined as Ω = (ωh ρ/C44)/π, Ψ = kh/π, respectively. The resulting dispersion curves are given in Figure 2. It clearly shows that the numerical results obtained by the present polynomial approach agree well with those obtained by the spectral collocation method, which validates our approach and program. The material in the above verification example is isotropic. We calculate the dispersion curves of Lamb wave in an orthotropic plate and make a comparison with the available results in literature [28] from the reverberation-ray matrix method, which serves as a further validation of our approach. The material is PZT-4, and the material parameters are listed in Table 2. Figure 3 shows the obtained frequency spectra. Here again, the agreement is quite good between our results and those from the reverberation ray matrix method. The material in the above verification example is isotropic. We calculate the dispersion curves of Lamb wave in an orthotropic plate and make a comparison with the available results in literature [28] from the reverberation-ray matrix method, which serves as a further validation of our approach. The material is PZT-4, and the material parameters are listed in Table 2. Figure 3 shows the obtained frequency spectra. Here again, the agreement is quite good between our results and those from the reverberation ray matrix method.  Then we discuss the convergence of the present polynomial approach. We present dispersion curves of propagating Lamb-like wave in a linear FGPPM plate with electric and magnetic open circuit and h = 1 mm, when the truncation order M takes 7, 8, 9 and 15, respectively, as shown in Figure 4. It can be seen that more and more order modes converge as M increases. When M = 7, the first three modes are convergent. The first four when M = 8, and the first seven when M = 9. So, we can think that at least the first (M − 1)/2 modes are convergent. Similarly, this can be concluded for the purely imaginary modes, and we don't present the dispersion curves of purely imaginary branches for saving space. For evanescent Lamb-like waves, we tabulate the results in Table 3 since graph is not convenient for comparison. These numerical results also show that the complex solutions are convergent as M increases. When M = 10 and M = 11, the first three modes are convergent. The first four when M = 12, the first five when M = 13, and the first six when M = 14. Obviously, the real solution is easier to converge than the complex one. From these results, good convergence of the present approach can be observed. We take M = 30 in this paper. Then we discuss the convergence of the present polynomial approach. We present dispersion curves of propagating Lamb-like wave in a linear FGPPM plate with electric and magnetic open circuit and h = 1 mm, when the truncation order M takes 7, 8, 9 and 15, respectively, as shown in Figure 4. It can be seen that more and more order modes converge as M increases. When M = 7, the first three modes are convergent. The first four when M = 8, and the first seven when M = 9. So, we can think that at least the first (M − 1)/2 modes are convergent. Similarly, this can be concluded for the purely imaginary modes, and we don't present the dispersion curves of purely imaginary branches for saving space. For evanescent Lamb-like waves, we tabulate the results in Table 3 since graph is not convenient for comparison. These numerical results also show that the complex solutions are convergent as M increases. When M = 10 and M = 11, the first three modes are convergent. The first four when M = 12, the first five when M = 13, and the first six when M = 14. Obviously, the real solution is easier to converge than the complex one. From these results, good convergence of the present approach can be observed. We take M = 30 in this paper. branches for saving space. For evanescent Lamb-like waves, we tabulate the results in Table 3 since graph is not convenient for comparison. These numerical results also show that the complex solutions are convergent as M increases. When M = 10 and M = 11, the first three modes are convergent. The first four when M = 12, the first five when M = 13, and the first six when M = 14. Obviously, the real solution is easier to converge than the complex one. From these results, good convergence of the present approach can be observed. We take M = 30 in this paper.

Full Dispersion Curves of Lamb Wave
Propagating waves have received a lot of attention, and here we put the emphasis on evanescent waves. We plot the full dispersion curves in 3D frequency-complex wave number space for a clearer visualization of the solutions and a better understanding of the nature of the modes, when necessary, with a different color for clarity. Figure 5a plots the full dispersion curves of Lamb wave for a linear FGPPM plate with electric and magnetic open circuit. Since the eigenvalues are computed for one ω at a time, the dispersion curves are constructed of unconnected dots and the points near the cut-off frequencies become sparse. We can observe that purely real and purely imaginary solutions appear in pairs of opposite signs and the complex ones appear in quadruples of complex conjugates and opposite signs. Purely real wave numbers correspond to the propagating wave, and purely imaginary and complex wave numbers correspond to the evanescent wave. For a given frequency, a certain small number of real branches exist together with an infinite number of complex and purely imaginary branches (mostly imaginary with few complexes in the given range). For clarity, Figure 5b shows one quadrant dispersion curves in a small range. For complex branches, most of them start from 0 frequency and end at the minima of the purely real branches. Occasionally, one connecting two purely imaginary branches appears. The real part of the complex branches is usually small. For purely imaginary branches, most of them start from 0 frequency and end at cut-off frequencies with increasing frequency, and some with small wave numbers start from one cut-off frequency and terminate the other one.
propagating wave, and purely imaginary and complex wave numbers correspond to the evanescent wave. For a given frequency, a certain small number of real branches exist together with an infinite number of complex and purely imaginary branches (mostly imaginary with few complexes in the given range). For clarity, Figure 5b shows one quadrant dispersion curves in a small range. For complex branches, most of them start from 0 frequency and end at the minima of the purely real branches. Occasionally, one connecting two purely imaginary branches appears. The real part of the complex branches is usually small. For purely imaginary branches, most of them start from 0 frequency and end at cut-off frequencies with increasing frequency, and some with small wave numbers start from one cut-off frequency and terminate the other one.   Figure 6 shows the phase velocity dispersion and attenuation curves of the first three propagating and complex branches. The dimensionless phase velocity and frequency and attenuation are defined by Vp = ω/(Re(k) · C 55 /ρ), fh = ωh/(2π C 55 /ρ) and Im(kh). We can find from these curves that the phase velocity of a propagating mode is decreased and gradually tends to a steady value with increasing frequency, but the velocity of an evanescent mode becomes bigger as well as the attenuation decreases. At high frequency, the evanescent mode has a very small attenuation, and its phase velocity is noticeably bigger than that of a propagating mode. For example, at fh = 2-3, the phase velocity of the second evanescent mode is about 8, but that of the propagating mode is below 2. Also, the wave dispersion is quite weak in this frequency range.
and Im(kh) . We can find from these curves that the phase velocity of a propagating mode is decreased and gradually tends to a steady value with increasing frequency, but the velocity of an evanescent mode becomes bigger as well as the attenuation decreases. At high frequency, the evanescent mode has a very small attenuation, and its phase velocity is noticeably bigger than that of a propagating mode. For example, at fh = 2-3, the phase velocity of the second evanescent mode is about 8, but that of the propagating mode is below 2. Also, the wave dispersion is quite weak in this frequency range.

Influences of Graded Field on Dispersion Curves and the Electromechanical Coupling Factor
Considering two graded shapes, cubic and sinusoidal graded fields. Figure 7 shows their dispersion curves of Lamb wave. The results show that the effect of the graded field on dispersion characteristics of Lamb wave is significant, including the propagating modes and evanescent modes. Comparison between Figure 7 and Figure 5b, we can notice that the imaginary part of the

Influences of Graded Field on Dispersion Curves and the Electromechanical Coupling Factor
Considering two graded shapes, cubic and sinusoidal graded fields. Figure 7 shows their dispersion curves of Lamb wave. The results show that the effect of the graded field on dispersion characteristics of Lamb wave is significant, including the propagating modes and evanescent modes. Comparison between Figures 7 and 5b, we can notice that the imaginary part of the complex branches for the sinusoidal graded case, at Ω = 0 plane, is bigger than that for the linear and cubic cases. Interestingly, for the cubic graded cases, the complex branch connecting two purely imaginary branches disappears and turns into a different one connecting a purely imaginary branch and a real branch.
For clarity, Figure 8 shows the frequency spectra and phase velocity spectra of Lamb propagating wave for the three graded fields. Obviously, the effect of the graded field is little on the low mode, but becomes significant with increasing the mode order and wave number. The phase velocity for the sinusoidal graded field is bigger than that for the linear graded field, while the linear bigger than the cubic. The reason lies in that the different graded fields result in different material volume distributions, and the wave velocity depends on the material properties. Figure 9 gives the variation curves of the three gradient fields in the z direction. The CoFe 2 O 4 content for the sinusoidal graded field is the highest, and the wave velocity of Ba 2 TiO 3 is slower than that of CoFe 2 O 4 .     The magneto-electromechanical coupling factor K 2 is an important parameter for designing acoustic wave devices. A high magneto-electromechanical coupling factor is expected in engineering applications. It is defined as [29]  The magneto-electromechanical coupling factor K 2 is an important parameter for designing acoustic wave devices. A high magneto-electromechanical coupling factor is expected in engineering applications. It is defined as [29] where Voc and Vsc are the phase velocities for the electric and magnetic open circuit and short circuit, respectively.
To illustrate the effect of graded field on the K 2 , we calculate the K 2 for S0 modes of four different FGPPM plates, as shown in Figure 10. We can find that the K 2 reaches a maximum at a certain wave number and tend to the same little value with increasing wave number, which implies the influence of the graded field on the energy propagation of Lamb wave in high-frequency zone is insignificant. It reaches a maximum from 4.4% for the sinusoidal graded field to 9.5% for the cubic graded field. They are located near kh = 2 and kh = 1.5 respectively. The K 2 for the cubic graded field is always bigger than that of the other three graded cases. Also the maximum of K 2 shifts to the smaller wave number when the graded power exponent is increasing. where Voc and Vsc are the phase velocities for the electric and magnetic open circuit and short circuit, respectively.
To illustrate the effect of graded field on the K 2 , we calculate the K 2 for S0 modes of four different FGPPM plates, as shown in Figure 10. We can find that the K 2 reaches a maximum at a certain wave number and tend to the same little value with increasing wave number, which implies the influence of the graded field on the energy propagation of Lamb wave in high-frequency zone is insignificant. It reaches a maximum from 4.4% for the sinusoidal graded field to 9.5% for the cubic graded field. They are located near kh = 2 and kh = 1.5 respectively. The K 2 for the cubic graded field is always bigger than that of the other three graded cases. Also the maximum of K 2 shifts to the smaller wave number when the graded power exponent is increasing. Figure 10. Effect of graded field on magneto-electromechanical coupling factor for S0 mode.

Wave Structure Analysis
The distributions of displacement and electric potential and magnetic potential fields can be obtained according to Equations (4) and (8). Considering a special position where the complex branch firstly collapses onto the real branch at about Ω = 1.0, as marked with a circle in Figure 5b. Figures 11 and 12 present the distributions of the physical quantities in the z and x directions when Ω = 1.01115, Ψ = 0.23216 − 0.04612i, and Ω = 1.01911, Ψ = 0.17938, respectively. As seen in these figures, the real branch propagates without any attenuation, and the complex branch exhibits an oscillatory distribution and propagates a very long distance, about a few tens of thicknesses of the plate. The displacement uz and electric potential and magnetic potential distributions change along the z direction in a nearly anti-symmetric manner. The displacement ux exhibits a nearly symmetric

Wave Structure Analysis
The distributions of displacement and electric potential and magnetic potential fields can be obtained according to Equations (4) and (8). Considering a special position where the complex branch firstly collapses onto the real branch at about Ω = 1.0, as marked with a circle in Figure 5b. Figures 11 and 12 present the distributions of the physical quantities in the z and x directions when Ω = 1.01115, Ψ = 0.23216 − 0.04612i, and Ω = 1.01911, Ψ = 0.17938, respectively. As seen in these figures, the real branch propagates without any attenuation, and the complex branch exhibits an oscillatory distribution and propagates a very long distance, about a few tens of thicknesses of the plate. The displacement u z and electric potential and magnetic potential distributions change along the z direction in a nearly anti-symmetric manner. The displacement u x exhibits a nearly symmetric manner. The distribution of displacement u z of the complex branch is very similar to that of the real branch, implying the evanescent wave mode converts into the propagating wave mode. Figure 10. Effect of graded field on magneto-electromechanical coupling factor for S0 mode.

Wave Structure Analysis
The distributions of displacement and electric potential and magnetic potential fields can be obtained according to Equations (4) and (8). Considering a special position where the complex branch firstly collapses onto the real branch at about Ω = 1.0, as marked with a circle in Figure 5b. Figures 11 and 12 present the distributions of the physical quantities in the z and x directions when Ω = 1.01115, Ψ = 0.23216 − 0.04612i, and Ω = 1.01911, Ψ = 0.17938, respectively. As seen in these figures, the real branch propagates without any attenuation, and the complex branch exhibits an oscillatory distribution and propagates a very long distance, about a few tens of thicknesses of the plate. The displacement uz and electric potential and magnetic potential distributions change along the z direction in a nearly anti-symmetric manner. The displacement ux exhibits a nearly symmetric manner. The distribution of displacement uz of the complex branch is very similar to that of the real branch, implying the evanescent wave mode converts into the propagating wave mode.

Merits of the Presented Method
Based on the above calculation of wave propagation in a FGPPM plate, we can summarize the following advantages of the presented method, which makes the method attractive.
(1) The complex mathematical issue is reduced to solve an eigenvalue problem, which is capable of accurately determining all the real, imaginary and complex solutions of a transcendental dispersion equation. (2) The conventional approaches (root-finding routines or finite element simulations) require an iterative search procedure or a far greater coding effort, to find complex roots. The present method can avoid tedious iterative two-variable search and is simple to program. It needs to

Merits of the Presented Method
Based on the above calculation of wave propagation in a FGPPM plate, we can summarize the following advantages of the presented method, which makes the method attractive.
(1) The complex mathematical issue is reduced to solve an eigenvalue problem, which is capable of accurately determining all the real, imaginary and complex solutions of a transcendental dispersion equation. (2) The conventional approaches (root-finding routines or finite element simulations) require an iterative search procedure or a far greater coding effort, to find complex roots. The present

Merits of the Presented Method
Based on the above calculation of wave propagation in a FGPPM plate, we can summarize the following advantages of the presented method, which makes the method attractive.
(1) The complex mathematical issue is reduced to solve an eigenvalue problem, which is capable of accurately determining all the real, imaginary and complex solutions of a transcendental dispersion equation. (2) The conventional approaches (root-finding routines or finite element simulations) require an iterative search procedure or a far greater coding effort, to find complex roots. The present method can avoid tedious iterative two-variable search and is simple to program. It needs to take a larger polynomial order to obtain solutions of the higher modes, which will cause more computer memory and long time. (3) The method is easy to implement and can be extended to complex structures such as multilayered or curved structures.

Conclusions
This paper presents an analytic method based on polynomial expansions for the determination of the full dispersion spectrum of the guided waves in the FGPPM plate. The correctness of the present method is verified via numerical comparison with available reference results. For the first time, the complete 3D dispersion curves of Lamb wave in a FGPPM plate are illustrated in a wide frequency range. The characteristics of the Lamb waves including the propagating and evanescent modes in various FGPPM plates are investigated. The emphasis on evanescent waves makes this work relevant for applications in the nondestructive evaluation of material or structural properties. Based on the above numerical results, some interesting conclusions can be drawn: (1) Superior to the conventional methods that necessitate an iterative search procedure to solve the complex roots of a dispersion equation, the presented analytic method can transform the set of differential equations for the acoustic waves into an eigenvalue problem in the form AX = kX to find the complex solutions. (2) Complex branches of the Lamb wave usually collapse onto the extremum of the real branches.
They exhibit both local vibration and local propagation, and some can propagate a quite long distance (more than ten times of the plate thickness). They will turn into the propagating modes with increasing frequency. (3) Some evanescent modes have a noticeably higher phase velocity than the propagating modes.
The phase velocity of the low order evanescent modes is more than four times larger than that of the propagating modes. Also, the wave dispersion of the evanescent mode is quite weak in a certain frequency range. (4) The magneto-electromechanical coupling factor of the guided wave in a FGPPM plate may be improved by adjusting the graded field. The coupling factor reaches a maximum from 4.4% for the sinusoidal graded field to 9.5% for the cubic graded field. The maximum of the magneto-electromechanical coupling factor for the S0 mode shifts to lower frequencies with increasing the gradient index.
Author Contributions: X.Z. and J.Y. proposed the studied problem and the corresponding solving method; X.Z. and Z.L. conducted the theoretical derivation and the computation; X.Z. and J.Y. analyzed and discussed the results; X.Z. and J.Y. wrote the paper.

Conflicts of Interest:
The authors declare no conflict of interest.