The Complex Rayleigh Waves in a Functionally Graded Piezoelectric Half-Space: An Improvement of the Laguerre Polynomial Approach

The research on the propagation of surface waves has received considerable attention in order to improve the efficiency and natural life of the surface acoustic wave devices, but the investigation on complex Rayleigh waves in functionally graded piezoelectric material (FGPM) is quite limited. In this paper, an improved Laguerre orthogonal function technique is presented to solve the problem of the complex Rayleigh waves in an FGPM half-space, which can obtain not only the solution of purely real values but also that of purely imaginary and complex values. The three-dimensional dispersion curves are generated in complex space to explore the influence of the gradient coefficients. The displacement amplitude distributions are plotted to investigate the conversion process from complex wave mode to propagating wave mode. Finally, the curves of phase velocity to the ratio of wave loss decrements are illustrated, which offers extra convenience for finding the high phase velocity points where the complex wave loss is near zero.


Introduction
Over the past decades, considerable attention has been focused on the investigation of functionally graded materials (FGMs), because the properties of this kind of heterogeneous composite materials can vary in preferred direction. These materials have been widely used in the thermal-protection systems of aerospace structures, so FGMs attract more and more interest of researchers in materials engineering. In recent years, however, FGPM has received rapt attention, because the material can be manufactured properly [1] with the development of material technology, and has been used in the manufacture of surface acoustic wave devices. For example, J. Qiu et al. [2] successfully used this kind of material to reduce the mechanical damping and heat generation of an actuator in 2003. All these cases show that FGPM is a promising material. On the other hand, surface acoustic wave technique has been adopted normally in signal processing, geophysics, and flow measurement, which makes it ever-increasingly important to study the wave propagation behaviors and characteristics in FGPM structures.
Many ground-breaking studies have been carried out by many researchers to explore the wave propagation behavior in FGPM, and various numerical methods have been proposed at the same time [3]. Early in 2004, Li et al. [4] analyzed the dispersion and anti-dispersion of Love waves in an FGPM half-space using the Wentzel-Kramers-Brillouin (WKB) method. Later, the WKB method was employed to investigate the propagation of Rayleigh surface waves in transversely isotropic FGPM half-space [5], and this method was also applied to study the dispersion relations and the stress distributions in an FGM layered half-space with initial stress [6]. Another approach to deal with the wave propagation in FGPM is the analytical technique. Du et al. [7] obtained the dispersion relations of the Love waves in FGPM layered half-space using an analytical approach, assuming that all material properties have the same exponential function distribution. Qian et al. [8] investigated the effect of gradient coefficient on dispersion relation and phase velocity by studying the characters of transverse surface waves using an analytical technique. The Love waves propagation in the FGPM buffer layer between the upper layer and the substrate was studied using the power series technique [9], and then the method was deployed to investigate the propagation behavior of Lamb waves in functionally graded piezoelectric-piezomagnetic material (FGPPM) plate [10]. Salah et al. combined the ordinary differential equation method with the general stiffness matrix method to investigate the propagation of Love and Rayleigh waves in FGPM half-space (see [11,12]). Very recently, Kong [13] and Mondal et al. [14] studied the propagating characteristics of SH waves in FGPM layered half-space using different method, respectively. On the problem of wave propagation in FGPM plates, C.Othmani et al. proposed a Legendre polynomial method to study the effect of the gradient coefficients on the Lamb phase velocity [15]. Later, they investigated the relationships of volume fraction, gradient coefficient, dispersion curves and phase velocity in a more complex structure of GaAs-FGPM-AlAs sandwich plate [16].
The above-mentioned studies on FGPM were conducted for the propagating waves. However, the studies on complex Rayleigh waves in the FGPM structures are rarely ever found. Zhang et al. [17] mentioned the complex acoustic parts of the Rayleigh wave in a half-space, and Rose [18] gave the explanation of the concept of complex waves, but their further research about the complex Rayleigh waves in FGPM half-spaces cannot be found. As is known, the wave equation can be expressed by angular frequency and wavenumber. When the equation is solved in complex space, it will generate three kinds of values: purely real value, purely imaginary values, or complex values. The real solutions represent propagating modes which do not need explanation anymore, and the solutions of purely imaginary values correspond to non-propagating modes because their amplitude decays exponentially. Nevertheless, it is worth noting that the non-propagating modes, corresponding to the complex values, exhibit the decay following an exponentially damped trigonometric distribution with distance [19]. There exist two big reasons the complex wave modes have been focused on: (1) these complex wave modes will travel a long distance if their imaginary parts are small enough, and they even could be converted to propagating wave modes when some conditions are satisfied; and (2) the complex wave mode can present much higher phase velocity than the surface acoustic wave [20], and further can be introduced to fabricate the SAW devices operating at higher frequencies.
Even though studies on complex Rayleigh waves in FGPM half-space are hardly found, the work on complex propagating waves has been carried out in other structures. It has been proved that some complex Lamb-like waves can be converted to propagating waves [21,22], but no research has been performed regarding the complex Rayleigh waves in FGPM half-space. In 1969, Lim and Farnell [23] began exploring the existence of the wave mode with phase velocity much higher than bulk wave. From then on, it became an ongoing important topic to study this kind of high phase velocity existing in the form of complex wave modes. Boyer et al. [24] solved the difficulty of poor estimation of the electroacoustical admittance for complex surface acoustic wave transducers by using fitting analytical function method. Benetti et al. not only carried out the theoretical investigation of complex surface acoustic wave properties in a layered structure [25], but also verified their results using experiment [26]. More research in this field can be found in the literature [27][28][29][30]. Thus far, most studies of complex surface acoustic waves have been limited in the structure of composite plate or simple layered half-space, and have not involved the structure of FGPM half-space.
Thus, in the present contribution, an improved method based on the Laguerre polynomial technique is proposed to investigate the complex Rayleigh waves in the FGPM half-space. The Laguerre polynomials are widely used to solve wave propagation problems in half-space because they are orthogonal over [0, ∞). Datta and Hunsinger [31] achieved surface wave solutions in piezoelectric half-space by using the Laguerre orthogonal polynomials while avoiding iterative calculation. [32] investigated the calculation of wave properties in an inhomogeneous half-space using the Laguerre orthogonal polynomials. Kim and Hunt investigated the surface acoustic wave propagation in multilyered half-space by using the Laguerre polynomial technique [33] and extended the method to examine the acoustic fields and velocities for surface acoustic wave propagation [34]. The conventional Laguerre orthogonal polynomial approach takes the wave numbers k as independent variables while taking the frequencies ω 2 as eigenvalues, but only one square term is related to ω in the Rayleigh wave governing equation. Thus, it can only obtain the real wavenumber solutions of propagating waves. However, the improved method in this paper takes the wave number k as eigenvalues while taking ω as independent variables. Both the linear and square terms of k exist in its governing equation, so the improved method not only can present the real wave number solution but also can get the solution of purely imaginary and complex wave number. The three-dimensional dispersion curves in electrically open cases are provided and discussed in order to analyze the influence of gradient coefficient. The distributions of displacement amplitudes are plotted and further applied to investigate the characteristics of complex Rayleigh waves. Finally, the curves of phase speed to the ratio of wave loss decrements are drawn, which provide extra convenience for finding the highest phase velocity points where the complex wave loss near zero. The results of this study can be used as the guidance of experimental measurement of Rayleigh surface waves in FGPM half-space.

Statement of the Problem
Consider an orthotropic FGPM elastic half-space which contacts with vacuum, and Rayleigh surface wave propagate along the x-axis as shown in Figure 1. The horizontal x-y coordinate plane is placed on the top surface of the media, and the z-axis of the Cartesian coordinate system is normal to the x-y plane; it is also the poling direction of the piezoelectricity. For the harmonic wave propagation considered in present paper, the body forces are assumed to be zero; thus, according to the linearity assumption, the dynamic equation for the half-space can be expressed as: where the comma subscript denotes a spatial coordinate, the double dot over u i denotes the second derivative with respect to time. The convention that j ≡ ∂ ∂x j is being used, and summation over a repeated subscript is understood. ρ is the mass density of the medium, u i is the mechanical displacement components in the ith direction, and σ and D are the mechanical stress tensor and electrical displacement tensor, respectively, which are related to the strain ε kl and the electric field intensity components E k in the constitutive relations as follows: where C ijkl is the elastic constants of physical medium, e kij is piezoelectric constants, and jk is dielectric constants. Their values vary continuously along the depth direction because the material is made up of FGPM, namely the material properties are the functions of depth in the z-axis. According to elasticity theory, the Cauchy's strain tensor is ε = 1 2 (u∇ + ∇u), and its six components are as follows: Based on the quasi-static Maxwell equation, the relationship of the electric field intensity E i and the electrical potential φ can be expressed as: In the present paper, only the coupled Rayleigh waves propagating along x-axis in the FGPM half-space are considered, and the poling direction is in z-direction, thus the mechanical and electrical boundary conditions (electric open case) at the free surface can be written as: Now, the problem of Rayleigh wave propagation in the FGPM half-space becomes that of finding the solution of the dynamic Equation (1) under the boundary condition of Equation (5).

Solution of the Problem
For the FGPM half-space polarized in the z-direction, the Rayleigh waves propagate along x-axis, so all field variables are independent on coordinate y. Thus, the governing Equation (1) in the sagittal plane are given by: To satisfy the boundary condition mentioned above, the unit step function is introduced here: then the mechanical and electric boundary conditions described in Equation (5) can be automatically incorporated into the constitutive equations of the half-space.
For the FGPM half-space in the present paper, whose material properties vary along the positive direction of z-axis, the elastic and other material coefficients can be written as: to obtain a more universal equation. Here, g 1 (z), g 2 (z), g 3 (z), and g 4 (z) are assumed to different gradient functions acting on the mass density ρ, elastic constants of physical medium C, piezoelectric constants e, and dielectric constants , respectively.
For a free harmonic half-space surface wave propagation in the x-direction, the displacement components and the electric potential can be expressed as: where U(z) and W(z) represent the amplitudes of vibration corresponding to the directions of x and z; X(z) represents the amplitude of electric potential; i is the imaginary unit;, ω is the angular frequency; and k is the complex wave vector in the propagation direction.
Equations (8)-(10a-c) are inserted into Equation (6) before dealing with them as the follows: Firstly, swap the two sides terms of the equations for the convenience of expression. Then, reorganize long expressions of every equation based on k's order and pick all the terms of same order with respect to k together. Finally, the governing differential equations in terms of displacement components become: In the above equations, the superscripts ( ) and ( ) are, respectively, the first and second partial derivatives with respect to z. U, W, and X correspond to the U(z), W(z), and X(z) in Equation (10a-c).
To get the solution of Equation (11a-c), U(z), W(z), and X(z) are expanded here into the products of Laguerre orthogonal polynomial series as follows: where p i m (i = 1, 2, 3) represent the expansion coefficients which correspond to displacement distributions for wave modes, L m (z) = e − z 2 P(m, z)/m!, and P(m, z) is the mth Laguerre polynomial, thus the function L m (z) forms a complete orthonormal set in the range (0, ∞), and is therefore suitable for calculating in half-space media. Even though the variable z is defined over (0, ∞) , the order of the expansion is always truncated to some value m = M in practice because the solution will converge within a finite number of terms.
Firstly, insert Equation (12a-c) into Equation (11a-c). Then, multiplying both sides of the four equations obtained in the last step by the complex conjugate L * j (z) with j running from 0 to M, 3 × (M + 1) equations appear. Finally, integrating the 3 × (M + 1) equations over z from 0 to infinity, and reordering these equations based on k 2 , k 1 and k 0 , achieves their matrix representation as follows: where P is a vector composed of p 1 , p 2 , and If only the wave number k is considered to be purely real for propagating waves, the solution can be obtained with ease through the conventional approach of the previous study [31]. The conventional method simplifies Equation (13) as the derivation of Equation (14a-c) in three steps, and Equation (14c) tends to be their final expression, which is a characteristic equation with ω 2 as the eigenvalues. Obviously, the method only can resolve the surface wave propagating problem in the real number range, because it takes the wavenumber as the independent variables and the angular frequencies only exist in real number.
However, in the present paper, the wave number k is non-real for complex wave modes, so the process to get the solution becomes more complicated because it entails calculating the propagating and non-propagating wave modes at the same time. Another difficulty lies in the fact that there exists a square term for k in the matrix in Equation (13), which will cause failure to obtain solutions through conventional eigenvalue methods. Thus, it is necessary to reform Equation (13) to become a classical eigenvalue problem in which the eigenvalue k appears linearly.
A new column vector is introduced here: which can be simplified into Substituting Equation (15) into Equation (13), the result becomes: After reorganizing, the new form of Equation (16) is Combining Equations (15) and (17) derives where E is an identity matrix which has the same dimension with A −1 B, R is a column vector composed of the p i m and q i m from Equation (15), The eigenvalues of coefficient matrix of left R in Equation (18) yield the complex wave number k, thereby the phase velocity ω/k, and the eigenvectors allow the calculation of the field distributions. Thus, solving the problem of complex Rayleigh wave propagation becomes calculating the eigenvalue problem of a specified matrix.

Numerical Results
To demonstrate the influences of the gradient coefficients on the characteristics of surface waves propagating on the free surface of FGPM half-space, some numerical analysis results are given in this section based on the solution derived in the previous section. Even though all the gradient functions g i (z) may be different because they are bound with different material properties, here, for the convenience of calculation, it is assumed that the variations of the elastic, piezoelectric, and dielectric coefficients and the mass density are the same along the z-axis direction, and the FGPM properties are controlled by the same gradient function g(z) = (1 + z) n , where n is the gradient coefficient. Therefore, the material coefficients can be expressed as: Lead zirconate titanate (PZT) is a commercially important piezoelectric material, and numerous studies on PZT-based piezoceramics have been conducted by researchers, thus, in the present paper, the set of FGPM half-space is assumed PZT-4, and the material properties (at z = 0) for the FGPM half-space are listed in the Table 1. Table 1. Properties of the piezoelectric PZT-4 [35] at z = 0.

Approach Validation
To verify the correctness of foregoing formulations, first the displacement curves using the present approach in a non-graded PZT-4 half-space are given, and then a comparison between the results and the published results from the literature [36] is provided. All these displacement curves are normalized with respect to w at z = 0, and given in Figure 2. Figure 2a was obtained from the present method and Figure 2b was taken from the literature. As can be seen from the comparison, the results from the method proposed in this paper and that from the literature agree very well.

Dispersion Curves for the FGPM Half-Space
For the convenience of understanding the nature of all wave modes and the straightforward visualization of solutions, the dispersion curves are presented in 3D graphics with different colors for clarity. When an angular frequency ω is given based on Equation (18), a vector is generated which consists of a series of wavenumbers including purely real, purely imaginary, and complex values. After constructing a three-dimensional coordinate system using frequency, purely real values, and purely imaginary values, the dispersion curves are obtained.
The first case is that the gradient coefficient n is 0.5 in the gradient function, and the resulting three-dimensional frequency spectra is illustrated in Figure 3. The graphic shows that the pure real (blue dotted curves in Figure 3), imaginary (red dotted curves in Figure 3), and complex roots (green dotted curves in Figure 3) always appear in pairs and are symmetrical about the coordinate plane. This indicates that, if a or ai is a solution, then −a or −ai is also a solution, and all complex solutions are evenly distributed in the four quadrants above in the form of conjugates and opposite signs. It is can be found that the blue curves from the pure real roots represent propagating wave modes, the red curves from the imaginary roots and the green curves from the complex roots correspond to the non-propagating wave modes. For pure real branches, the first branch starts from the cut-off frequency of the first imaginary branch, which is very different from Lamb-like waves starting at zero, but the other real branches' tendency are similar to Lamb-like waves. As for pure imaginary branches, the first three branches are different from the others because the first branch starts at the frequency zero while ending at the starting point of the real branch, and their next two branches exist between two adjacent cut-off frequencies. Except for the above, the other imaginary branches after the third branch start at a certain cut-off frequency and linked in three-dimensional space by the green complex branches, some of them are coupled together very well. For the green complex branches, they start at frequency 0 and end at the minimum points of the red pure imaginary or the blue pure real branches. It is worth mentioning that the first two complex branches are non-propagating modes at small frequencies and decay more slowly as frequency f increase, but their endpoints falls on the blue complex branches. This means the complex branches degenerate into pure real branches, i.e. these non-propagating modes are converted into propagating modes at these points as the complex curves proceed.
To investigate the influence of the heterogeneous character of the FGPM on the dispersion curves, the gradient coefficient n is chosen to be 1 and 2 in the present numerical analysis, and a dispersion comparison with n is 0.5 is presented. As mentioned above, all these curves are symmetrical about the coordinate plane, thus only the dispersion curves in of one quadrant are given. Figure 4, with the same color scheme as in Figure 3, illustrates two three-dimensional dispersion curves for Rayleigh waves in FGPM half-space with different gradient coefficients. The results in Figures 3 and 4 show that the frequency of the same order propagating mode increases as the gradient coefficient n increases for a specified wavenumber. This means that the speed of Rayleigh wave propagation becomes faster as the gradient coefficient n increases because the propagation speed of any point on the blue real branches is proportional to the ratio of f /k. For pure imaginary modes, the high order branches tend to become more independent of each other as the gradient coefficient n increases, and wave mode conversions happen between them and complex waves. However, the most interesting thing happens with complex waves. Except for those local short complex modes appearing between high order pure imagine modes, the endpoints of complex modes fall on the pure real branches and move forward the f -axis gradually as the gradient coefficient n increases. As a result, this kind endpoints is reduced by one when n is 2, but the other falls on the pure imaginary branch, as illustrated in the Figure 4.

Displacement Amplitude Distributions for the FGPM Half-Space
Next, the displacement amplitude distributions is given to investigate the differences between the propagating waves and the complex Rayleigh waves. Taking the expression k = k R + ik I (k R ∈ R, k I ∈ R) as the conventional solution of the wave number, substituting the expression into Equation (12a-c), and then inserting the result into Equation (10a-c), the new form of displacement can be written as: Considering Euler's formula and assuming t = 0, the displacement can be calculated by the following, The displacement field distributions are represented in Figure 5 evaluated at f = 22.3825 in the thickness direction (z) and the propagation direction (x), respectively. Figure 5a shows that the damped exponential distribution happens to the pure imaginary wave modes corresponding to those imaginary wave numbers. The complex Rayleigh wave modes exhibit damped sinusoidal filed distribution with rapid decay, as shown in Figure 5b. This phenomenon becomes more pronounced when the ratio of the real part to the imaginary part in the complex solution is relatively large, as demonstrated in Figure 5c. As the ratio of the real part to the imaginary part becomes larger, the complex wave mode must be converted into pure real wave mode corresponding to propagation wave, as demonstrated in Figure 5d. This consequence is verified through some complex wave branches, such as the curve CA in Figure 4.

The Effects of Gradient Coefficients on the High Phase Velocity Points
To find the highest phase velocity points where the imaginary part of complex wave numbers are near zero is one purpose of the present study. Furthermore, the effects of gradient coefficients on these points is examined. Thus, it is essential to investigate the phase velocity curves of complex wave numbers with different gradient coefficient in the FGPM half-space. Firstly, the real part and imaginary part of complex wave numbers are set to horizontal coordinates that are perpendicular to each other. Then, the phase velocity is set to the vertical coordinates. As a result, a three-dimensional coordinate system about phase velocity and complex wave number is established. The phase velocity curves of complex wave numbers with different gradient coefficient are illustrated in Figure 6, in which the purely imaginary branches are removed to clearly demonstrate the conversion process from complex wave modes to real wave modes. It can be seen clearly that some complex wave modes approach purely real propagating wave modes in this figure, and the phase velocity of the first point which imaginary part near zero becomes large with the increase of gradient coefficient.
To find the exact points where the complex wave modes are near purely real propagating wave modes, the loss ratio η is given. Here, denote the ratio η = Imaginary(k)/Real(k), and set it to the abscissa. Then, set the phase velocity to ordinate, so that the coordinate system is established about loss ratio and phase velocity, as illustrated in Figure 7. As can be seen, the phase velocity of the first conversion point is 6.6 km/s when the gradient coefficient n = 1. However, the velocity increases to 7.4 km/s when the gradient coefficient n = 2. Thus, it can be concluded that the changes of the gradient coefficients have a great effect on the first point where the imaginary parts of the complex wave modes approach zero, and where the phase velocities of complex waves travel much more quickly than the propagating surface waves. The calculated velocity of propagating surface wave in PZT-4 is 2.2 km/s [36] and the measured value is 2.12 km/s [37].

Conclusions
A polynomial expansion approach for solving the both propagating and non-propagating Rayleigh modes in the free surface of FGPM half-space has been developed, and the full dispersion spectrum and the displacement amplitude field distributions have been investigated by the present method while presenting all the possible solutions of the pure real, pure imaginary, and the complex. As a result, the following conclusions can be drawn:

1.
The improved method proposed in this paper can be used to solve the problem of complex Rayleigh waves in the FGPM half-space. This method can not only provide real roots, but also imaginary and complex roots.

2.
The curves of dispersion and displacement demonstrate that some complex wave branches could be converted to the propagation waves.

3.
The gradient coefficient n has great effect on the frequency of the same branch dispersion curve, and the corresponding frequency will increase with the n increase.