Asymptotic Solution and Numerical Simulation of Lamb Waves in Functionally Graded Viscoelastic Film

To investigate Lamb waves in thin films made of functionally graded viscoelastic material, we deduce the governing equation with respect to the displacement component and solve these partial differential equations with complex variable coefficients based on a power series method. To solve the transcendental equations in the form of a series with complex coefficients, we propose and optimize the minimum module approximation (MMA) method. The power series solution agrees well with the exact analytical solution when the material varies along its thickness following the same exponential function. When material parameters vary with thickness with the same function, the effect of the gradient properties on the wave velocity is limited and that on the wave structure is obvious. The influence of the gradient parameter on the dispersion property and the damping coefficient are discussed. The results should provide nondestructive evaluation for viscoelastic material and the MMA method is suggested for obtaining numerical results of the asymptotic solution for attenuated waves, including waves in viscoelastic structures, piezoelectric semiconductor structures, and so on.


Introduction
Lamb waves, which are a type of plain strain wave in a thin film or a plate with a traction-free boundary, are widely used in nondestructive evaluation. Early research reported on Lamb waves focused on isotropic elastic plates [1]. Since then, scientists have directed more attention to Lamb waves in plates made of various materials, including viscoelastic materials [2], functionally graded materials (FGMs) [3], piezoelectric materials [4], and piezoelectric-piezomagnetic materials [5]. To detect material properties or damage to the structures, much research has been focused on guided waves in composite structures based on numerical and experimental methods [6,7].
FGMs were proposed by scientists as a kind of thermal-protection material in the 1990s [8]. In FGM structures, the material parameters are not constant and vary along one direction continuously. With the development of material technology, the FGM technique has been used not only for common elastic material but also for some smart materials, including piezoelectric [9,10] and piezoelectric-piezomagnetic materials [11]. To evaluate the mechanical properties of FGM structures, researchers have investigated various elastic waves in FGM structures, such as Lamb waves, horizontal shear (SH) waves [12], Love waves [13], and Rayleigh waves [14].
To address wave propagation problems in heterogeneous media, both numerical and analytical methods are employed for solving the wave-governing differential equations with variable coefficients. The main idea of numerical methods is to divide the functionally graded material into multilayer models and to simplify each sublayer as a homogenous layer [15][16][17][18]. Scientists have also proposed some analytical solutions for wave propagation problems in different heterogeneous structures. These methods include exact analytical expressions for material parameters following the same exponential function [19], the Wentzel-Kramers-Brillouin (WKB) method for large-wave-number [20] or cutoff problems [21,22], and the special function method for material parameters following some special function [23]. In recent decades, researchers suggested that these equations can be solved by using a power series method [11,24] and a Legendre polynomial method [25,26], which are fit for solving the wave propagation problem in heterogeneous structures in arbitrary cases in which material parameters vary continuously and slowly. The form of the dispersion equations based on these two methods contains series items. Therefore, these dispersion equations should be solved numerically.
It is found that not only elastic materials but also viscoelastic materials in nature have gradient properties. For example, when a material undergoes subsurface aging or subsurface damage, the elastic modulus varies along the thickness of the damaged subsurface region and mechanical gradient characteristics appear [27]. This should also occur for viscoelastic materials. For the wave propagation problem in viscoelastic structures, Lu et al. [28] found that the attenuation of Lamb waves increases with the increase of the thickness of the viscoelastic layer and that the mode is transformed as well. Compared with the propagation characteristics of Love waves in an elastic medium, the energy of Love waves in a Kelvin-Voigt viscoelastic medium is obviously attenuated, as shown by Zhang et al. [29]. SH waves have one displacement component. Yu et al. [30] deduced the dispersion equations for SH waves in orthotropic viscoelastic hollow cylinders. There are few studies on the propagation of Lamb waves with two displacement components in viscoelastic complex structures, and most of them use the Legendre polynomial method [31].
The dispersion equations for wave propagation in a viscoelastic material comprise a set of complex coefficient transcendental equations. To solve the transcendental equations with complex variables, Qian et al. [32] comprehensively analyzed the applicability of the parabolic Newton iteration method, the binary dichotomy method, and the modulus value convergence method. However, when the power series method is employed to solve the wave propagation problem in a functionally graded viscoelastic material (FGVM) structure, the dispersion equation, which is a transcendental equation with complex numbers in series form, is difficult to solve based on the above numerical simulation method. For example, the Newton iteration method requires that the solution be in the form of a display function rather than a series, while the binary dichotomy method and the modulus value convergence method might lead to the existence of spurious solutions.
In this study, we investigate the dispersion and attenuation characteristics of Lamb wave propagation in a thin film made of FGVM, which follows the Kelvin-Voigt model [33]. The governing equations with a displacement function are deduced and the power series asymptotic solution is obtained by using the power series method. Because the series has no explicit expression for the function, we propose the minimum module approximation (MMA) method for solving the complex coefficient dispersion equation. The detailed process of the MMA method, the existence analysis of its solution, and its optimization are given. The reliability of the power series solution is verified by comparison with the exact analytical solution for Lamb wave propagation in a functionally graded viscoelastic film. The dispersion and attenuation characteristics of Lamb wave propagation under different gradient parameters are discussed, and the damping coefficients are analyzed. Conclusions based on these results can provide a theoretical basis for nonhomogeneous viscoelastic structure nondestructive testing.

Basic Equation for Lamb Waves in FGVM Film
Consider Lamb waves propagating in an isotropic functionally graded viscoelastic film along the x direction, as shown in Figure 1. The thickness of the film is h. The z direction is along the thickness direction. Let u, v, and w represent the displacement in the x, y, and z directions, respectively. For Lamb waves propagating is this structure, the displacement should satisfy:  If we let subscripts 1, 2, and 3 represent x, y, and z, respectively, then the stress tensor ij σ , In the Kelvin model, the constitutive equation can be expressed as: where , , , 1 3 i j k l = − , kl S is the strain tensor, η represents the viscosity coefficient, G is the shear modulus, and K is the bulk modulus. Equation (2) can be rewritten as: where ijkl c and ijkl c are components of the elastic tensor and viscosity tensor, respectively. In an FGVM, both these elastic parameters as well as the mass density ρ are not constants but are functions of z. The motion equations have the following form: If we let subscripts 1, 2, and 3 represent x, y, and z, respectively, then the stress tensor σ ij , i, j = 1 − 3, can be divided into two parts: the deviation stress σ ij and the spherical stress tensor σ m , where σ m = σ kk /3, k = 1 − 3, and the repeated index in the subscript implies summation with respect to that index.
In the Kelvin model, the constitutive equation can be expressed as: where i, j, k, l = 1 − 3, S kl is the strain tensor, η represents the viscosity coefficient, G is the shear modulus, and K is the bulk modulus. Equation (2) can be rewritten as: where c ijkl and c ijkl are components of the elastic tensor and viscosity tensor, respectively. In an FGVM, both these elastic parameters as well as the mass density ρ are not constants but are functions of z. The motion equations have the following form: The relation between the strain and the displacement deduced from Equation (1) is: By using index reduction, the original fourth-order elastic parameters can be rewritten to second-order elastic parameters. By substituting Equation (5) into the constitutive Equation (3), we obtain the following component forms of the stress: where the parameters in isotropic materials satisfy c 11 − c 13 = 2c 44 and c 11 − c 13 = 2c 44 . Substitution of Equation (6) into Equation (4) leads to the following governing equations with respect to the displacement components: Lamb waves propagating in a functionally graded viscoelastic film must satisfy not only the governing equations but also the traction-free conditions of the film, which are expressed as:

Power Series Solution
The solutions of the governing equations can be expressed as: where i is the imaginary unit, k and c are the wave number and wave velocity, respectively, and U(z) and W(z) are the unknown amplitudes of the displacement. Substitution of Equation (9) into Equation (7) leads to: whereĉ 44 = c 44 − iΩc 44 ,ĉ 11 = c 11 − iΩc 11 , Ω = ck is the frequency, and the prime symbol ( ) represents differentiation with respect to thickness z. Bothĉ 44 andĉ 11 are functions of z and Ω. Suppose that the material parameters vary along the thickness direction slowly, so that, for a certain Ω, the material parameters of the isotropic FGVM film can be expressed as follows: Suppose that the material parameters can be expressed in the power series form: Therefore, the solutions of Equation (10) can also be expressed in a power series form as: By substituting Equations (11)-(13) into Equation (10) and equating the coefficient of (z/h) n to zero, the following recursive equations can be obtained: n l=0 (l + 2)(l + 1)a 2 n−l s l+2 + n l=0 (n − l + 1)(l + 1)a 2 n−l+1 s l+1 + (kh) 2 n l=0 a 3 n−l c 2 − a 1 n−l s l where s 0 , s 1 , t 0 , and t 1 are undetermined coefficients. For l ≥ 2, all of the s l and t l are linear functions of s 0 , s 1 , t 0 , and t 1 .
To simplify calculating the relation between s l , t l and s 0 , s 1 , t 0 , t 1 , let: where j = 1-4 and I is a 4 × 4 unit matrix. Therefore, the solution of Equation (10) can be rewritten as: where the constants C j (j = 1-4) are to be determined. For n = 0 and 1, s nj and t nj are defined by expression (15); for other values of n, s nj , and t nj can be determined by solving Equation (14). By substituting (16) into the boundary conditions, we then obtain the following linear algebraic equations for determining constants C j (j = 1-4): The sufficient and necessary condition for the existence of a nontrivial solution is that the determinant of the coefficient matrix must vanish. Therefore, for the dispersion relation for Lamb waves, there exists: where where j = 1-4, and other items of T ij equal zero. The superscripts 0 and h represent the material parameters at the bottom and upper surfaces, respectively.
Owing to the existence of the complex relation, Equation (18) is a complex coefficient transcendental equation. In this paper, we suppose that k and the wavelength λ are both real numbers. The wave velocity c contains both real and imaginary parts. The real part of the wave velocity represents the phase velocity, and the imaginary part is related to the attenuation characteristic of the wave.

MMA Method
For solving an equation with a complex variable, we should obtain a solution for which both the real part and the imaginary part of the equation should be zero. Consider the complex equation: where x and y are the real and imaginary parts, respectively, of the complex variable z, and f, which is a function of z, is also complex. We suppose z = a + ib is the solution of Equation (19). It should satisfy the conditions: where Re[ f (x, y)] and Im[ f (x, y)] represent the real and imaginary parts of the function f (x, y). where Mod[ f (x, y)] is the module of the function f (x, y) and the square of the module is expressed by the function G(x, y).
Suppose that the solution of Equation (19) exists and is unique in the region (x 0 , x 0 ), y 0 , y 0 . The model in Figure 2 is used to illustrate the solution steps. The first loop step is shown in Figure 2a: We divide the solution region into n × n grids, obtaining (n + 1) 2 nodes in total, and calculate the module of each node by using Equation (21) and find the node (a 1, b 1 ) satisfying: where  represent the real and imaginary parts of the function where ( ) , f x y and the square of the module is expressed by the function G(x, y).
Suppose that the solution of Equation (19) exists and is unique in the region ( ) , y y . The model in Figure 2 is used to illustrate the solution steps. The first loop step is shown in Figure 2a: We divide the solution region into n × n grids, obtaining (n + 1) 2 nodes in total, and calculate the module of each node by using Equation (21) and find the node (a1, b1) satisfying: where For the second loop step, based on n, we obtain the four points (a1−Δx0,b1−y0), (a1+Δx0,b1−y0), (a1+Δx0,b1+y0), and (a1−Δx0,b1+y0), remesh the square determined by the four points as the vertices into n × n grids once again, and also calculate the module of each node to find the minimum.
For the Nth loop step (Figure 2b), we repeat the above step, obtaining: is the node with the minimum module after N time steps. The approximate solution of Equation (19) i N N z a b = + , can then be obtained.

Optimization of the MMA Method
The MMA method can be applied to solve equations with complex variables. To optimize the method, the following function q is introduced: For the second loop step, based on n, we obtain the four points (a 1 −∆x 0 ,b 1 −y 0 ), (a 1 +∆x 0 ,b 1 −y 0 ), (a 1 +∆x 0 ,b 1 +y 0 ), and (a 1 −∆x 0 ,b 1 +y 0 ), remesh the square determined by the four points as the vertices into n × n grids once again, and also calculate the module of each node to find the minimum.
For the Nth loop step (Figure 2b), we repeat the above step, obtaining: where (a N , b N ) is the node with the minimum module after N time steps. The approximate solution of Equation (19) z = a N + ib N , can then be obtained.

Optimization of the MMA Method
The MMA method can be applied to solve equations with complex variables. To optimize the method, the following function q is introduced: where q represents the average percentage of solution area reduction with each calculation and n 2 and (n + 1) 2 are the number of grids and nodes of the solution region, respectively.
To optimize the MMA method, we should find the n value needed to satisfy that q reaches its minimum. In other words, ln q should also be the minimum. Therefore, n should satisfy: The solution of Equation (25) is n = 3.77. Because the number of grids needs to satisfy the condition of positive integers, we take the approximate solution at n = 4.
In practical numerical analysis, we always predict with a certain level of uncertainly the region over which the equation with complex variables has a solution. If, in the first step, (a 1 , b 1 ) lies on the boundary of the region, (x 0 , x 0 ), y 0 , y 0 , the solution might not lie in the region. In this case, n should be selected to be a larger number to certify the existence of the solution. However, a spurious solution might exist if the calculation region is too large. Normally, the solutions of these problems are always irrational. This means that we can find the solution as the module infinites approaches zero. We should check for the convergence of solution by testing the ratio of the module reduction in several continuous steps. For example, we can calculate the ratio of the module for every three loops and judge the convergence to avoid a spurious solution.
It is worth noting that, if the minimum modulus is located at the boundary in the first calculation, there might be no solution in the computational domain. If the minimum modulus is not at the boundary after mesh refinement, the solution exists in the computational domain. Otherwise, the computational domain needs to be enlarged and recalculated. To verify the existence of the solution and avoid a spurious solution, we suggest that n should be selected as 6-8 in the first loop step in practical calculations.

Comparison with the Exact Analytical Results
To verify the validity of the power series method, the exact analytic solution and the asymptotic solution of the power series are compared when all material parameters vary with the same exponential function. The exact analytical solution can be obtained for waves propagating in the special FGVM thin film. The governing equations can be simplified to ordinary differential equations with constant coefficients, and the analytical solution can be obtained directly.
We suppose that the material parameters follow: where λ 0 , µ 0 , and ρ 0 are material parameters of the film lower surface at z = 0, η is the viscosity coefficient, ξ is a constant and is selected as 10 −5 , and p represents the gradient parameter. The analytical solution in this condition can be selected as a reference for the solution of the power series reported in this paper. The material parameters used in this paper can be deduced as: The displacement amplitude can also be expressed in an exponential function form as: By substitution of Equation (28) into Equation (10) the homogeneous linear equations for the undetermined coefficients A and B can be deduced as: Equations (29) comprise a set of linear homogeneous equations with respect to A and B. From the necessary and sufficient conditions for the existence of a nontrivial solution, we obtain that the determinant of the coefficient matrix is equal to zero: Considering that Equation (30) is a fourth-order equation, we suppose that the solution is α j (j =1-4). The relation between A and B is derived by calculating Equation (29) as follows: The displacement amplitude solution of Equation (10) can be rewritten as: Similarly, by considering the boundary conditions, we then obtain the dispersion equation: where In numerical analysis, the normalized wave velocityĉ and the dimensionless wave number kh are applied for describing the wave propagation property. The normalized dimensionless wave velocitŷ c satisfies:ĉ = c/c sh where c sh = G/ρ, which is the bulk shear wave velocity. The Poisson ratio is a constant and satisfies ν = 0.25. To evaluate the accuracy and precision of the power series and MMA methods, the relation between the normalized wave velocity and the dimensionless wave number for Lamb waves propagating in the special FGVM thin film is plotted in Figure 3. When p = 0, the FGVM thin film becomes a homogenous viscoelastic thin film. Figure 3a,b present the real and imaginary parts of the normalized wave velocity, respectively. It is found that the solution obtained by using the power series method agrees well with the exact analytical solution. By checking the results of the change of phase velocity, we find that there exists little difference between the wave velocity curves for Lamb waves in the homogenous viscoelastic thin film and in the special FGVM thin film. Normally, the dispersion curves of Lamb waves in homogenous film can be determined by bulk shear wave velocity and the Poisson ratio. Both the bulk shear wave velocity By checking the results of the change of phase velocity, we find that there exists little difference between the wave velocity curves for Lamb waves in the homogenous viscoelastic thin film and in the special FGVM thin film. Normally, the dispersion curves of Lamb waves in homogenous film can be determined by bulk shear wave velocity and the Poisson ratio. Both the bulk shear wave velocity and the Poisson ratio in homogenous thin film are same as those in the special FGVM thin film. It can be used to explain that the dispersion curves of the two cases are almost identical. This implies that we cannot measure the gradient parameters by variation of the wave velocity We further study the wave structure of Lamb wave propagation in different viscoelastic thin films. The normalized displacement amplitude is defined as: where U(0), which represents displacement component at z = 0, is selected to be 1 in the numerical analysis, Re is the real part of the complex number, and the sign function satisfies The normalized displacement amplitude of the first two modes at kh = π and kh = 2π are plotted in Figure 4. The curves obtained from the exact solution and these obtained by using the power series method coincide completely. In a homogenous viscoelastic thin film, the wave structure is symmetric or antisymmetric. However, because of the asymmetric properties of the FGVM thin film, the displacement amplitudes are not symmetric. By checking the results of the change of phase velocity, we find that there exists little difference between the wave velocity curves for Lamb waves in the homogenous viscoelastic thin film and in the special FGVM thin film. Normally, the dispersion curves of Lamb waves in homogenous film can be determined by bulk shear wave velocity and the Poisson ratio. Both the bulk shear wave velocity and the Poisson ratio in homogenous thin film are same as those in the special FGVM thin film. It can be used to explain that the dispersion curves of the two cases are almost identical. This implies that we cannot measure the gradient parameters by variation of the wave velocity We further study the wave structure of Lamb wave propagation in different viscoelastic thin films. The normalized displacement amplitude is defined as: , which represents displacement component at z=0, is selected to be 1 in the numerical analysis, Re is the real part of the complex number, and the sign function satisfies The normalized displacement amplitude of the first two modes at kh π = and 2 kh π = are plotted in Figure 4. The curves obtained from the exact solution and these obtained by using the power series method coincide completely. In a homogenous viscoelastic thin film, the wave structure is symmetric or antisymmetric. However, because of the asymmetric properties of the FGVM thin film, the displacement amplitudes are not symmetric. The relation between the gradient parameter and the ellipticity of particles at kh π = and 2 kh π = is plotted in Figure 5. It is found that the influence of the gradient property on the ellipticity of a To further investigate the influence of the gradient property on the displacement, we denote the ellipticity of particles on lower and upper surfaces as χ 0 = |w 0 /u 0 | and χ h = |w h /u h |, respectively. The relation between the gradient parameter and the ellipticity of particles at kh = π and kh = 2π is plotted in Figure 5. It is found that the influence of the gradient property on the ellipticity of a particle on the surface is more obvious than that on the wave velocity.  The relation between the gradient parameter and the ellipticity of particles at kh π = and 2 kh π = is plotted in Figure 5. It is found that the influence of the gradient property on the ellipticity of a particle on the surface is more obvious than that on the wave velocity.

Material Parameters Varying Linearly
For numerical analysis with the theoretical model described above, we assumed that the Lamé parameters λ and µ, mass density ρ, and viscosity coefficient η in the functionally graded viscoelastic film varied as follows: where λ 0 , µ 0 , and ρ 0 are material parameters of the film lower surface at z = 0; η is the viscosity coefficient; ξ is a constant and is selected as 10 −5 (except in Section 5.2.3); and p 1 , p 2 and p 3 are the gradient parameters of λ, µ, and ρ, respectively (0 ≤ p 1 , p 2 , p 3 < 1). The material parameters used in this paper can be deduced as:

All Material Parameters Varying Identically
Suppose that all material parameters vary along the thickness direction linearly and identically, i.e., p i = p(i = 1, 2, 3). The wave velocity plotted as a function of wave number when p = 0, 0.3, 0.5, and 0.7 is shown in Figure 6. By comparing the results for the case in Section 5.1, a similar conclusion can be reached. If all the material parameters vary along the thickness direction identically, the wave velocity curves almost coincide. 11 44

All Material Parameters Varying Identically
Suppose that all material parameters vary along the thickness direction linearly and identically, i.e., ( ) The wave velocity plotted as a function of wave number when p = 0, 0.3, 0.5, and 0.7 is shown in Figure 6. By comparing the results for the case in Section 5.1, a similar conclusion can be reached. If all the material parameters vary along the thickness direction identically, the wave velocity curves almost coincide. Similarly, the wave structures are also plotted in Figure 7. It is found that the gradient parameter has an obvious influence on the wave structure. It is also shown in Figure 7 that the ellipticity of particles on the upper and lower surfaces is different owing to the gradient property of the FGVM film. Considering testability, we suggest that the ellipticity of a particle on the surface can be applied for measuring the gradient parameter when all material parameters vary identically. In these cases, the bulk wave velocity, including the shear wave velocity c sh and the longitudinal wave velocity c L , where c L = (λ + 2µ)/ρ, are constants. This implies that, if the bulk wave velocities are constants, the wave velocity of Lamb waves in the FGVM thin film are almost similar to that in a homogenous film.
Similarly, the wave structures are also plotted in Figure 7. It is found that the gradient parameter has an obvious influence on the wave structure. It is also shown in Figure 7 that the ellipticity of particles on the upper and lower surfaces is different owing to the gradient property of the FGVM film. Considering testability, we suggest that the ellipticity of a particle on the surface can be applied for measuring the gradient parameter when all material parameters vary identically.

Material Parameters Varying Independently
To investigate the influence of the elastic modulus and density gradient on dispersion and attenuation characteristics of Lamb waves in gradient viscoelastic film, we chose three types of films

Material Parameters Varying Independently
To investigate the influence of the elastic modulus and density gradient on dispersion and attenuation characteristics of Lamb waves in gradient viscoelastic film, we chose three types of films for which the gradient parameters are: Film A is a homogeneous film, which can be used for referencing the propagation characteristics in the gradient film. The elastic modulus in film B is a constant, and the density increases along the thickness of the film and the density of the lower surface is the same as that of the homogeneous film. Conversely, the density in film C is a constant, and the elastic modulus varies along the thickness direction linearly.
The real and imaginary parts of the wave velocity in the three types of films are shown in Figure 8. When the mass density increases, the real part and the absolute value of the imaginary part of the wave velocity of each mode are less than these in the homogenous film; when the elastic modulus increases, the real part and the absolute value of the imaginary part of the dimensionless wave velocity of each mode are larger than these in the homogenous film.

Relative Viscosity Coefficient Varying Independently
In engineering application, the relative viscosity coefficient might vary along one direction because of the environment. However, in these cases, the mass density and the elastic parameters might not change. To reveal the influence of the gradient relative viscosity coefficient on the wave property, we suppose that the relative viscosity coefficient ξ varies along the thickness direction and that other parameters including λ , μ , and ρ are constants. The material parameters are: The wave velocity is plotted as function of wave number in Figure 9. In Figure 9a, the curves for the real part of the phase velocity almost coincide. This suggests that the influence of the gradient relative viscosity coefficient on the dispersion curves is too slight to measure. However, obvious differences can be observed in Figure 9b, which describes the relation between the imaginary parts of the wave velocity. The absolute value of the imaginary part of the dimensionless wave velocity increases with the increase of the gradient relative viscosity coefficient. The physical meaning of the imaginary parts of the wave velocity will be discussed in the next section.

Relative Viscosity Coefficient Varying Independently
In engineering application, the relative viscosity coefficient might vary along one direction because of the environment. However, in these cases, the mass density and the elastic parameters might not change. To reveal the influence of the gradient relative viscosity coefficient on the wave property, we suppose that the relative viscosity coefficient ξ varies along the thickness direction and that other parameters including λ, µ, and ρ are constants. The material parameters are: The wave velocity is plotted as function of wave number in Figure 9. In Figure 9a, the curves for the real part of the phase velocity almost coincide. This suggests that the influence of the gradient relative viscosity coefficient on the dispersion curves is too slight to measure. However, obvious differences can be observed in Figure 9b, which describes the relation between the imaginary parts of the wave velocity. The absolute value of the imaginary part of the dimensionless wave velocity increases with the increase of the gradient relative viscosity coefficient. The physical meaning of the imaginary parts of the wave velocity will be discussed in the next section. the real part of the phase velocity almost coincide. This suggests that the influence of the gradient relative viscosity coefficient on the dispersion curves is too slight to measure. However, obvious differences can be observed in Figure 9b, which describes the relation between the imaginary parts of the wave velocity. The absolute value of the imaginary part of the dimensionless wave velocity increases with the increase of the gradient relative viscosity coefficient. The physical meaning of the imaginary parts of the wave velocity will be discussed in the next section.

Influence of Gradient Parameter on Wave Attenuation
Equation (18) is a complex equation. As the wave number is a real number, the wavelength is also a real number, the obtained wave velocity is complex, and the product Ω of the wave velocity c and the wave number k is complex, which can be expressed as follows: where ω and ω are the real and imaginary parts of Ω, respectively, ω is frequency, and ω is related to the attenuation of the wave amplitude. From Equation (40), we have: In viscoelastic materials, the wave propagation process is essentially a quasi-periodic motion, and the period of the particle displacement is determined by the phase velocity. The period is expressed as follows: To analyze the attenuation trend, we define the amplitude ratio of the adjacent period as the damping coefficient γ, given by: In this study, the normalized product of frequency and thicknessωh is selected to be the abscissa. If The influence of the gradient properties on the damping coefficient is plotted in Figure 10. The damping coefficient increases with the increase of the frequency. When material parameters are constants, or material parameters vary along the thickness direction with the same exponential function, or both Lamé parameters and mass density vary linearly, the damping coefficient of Mode 1 and Mode 2 is similar, as shown in Figure 10a,b. However, when the Lamé parameters and mass density do not vary identically, a difference in the damping coefficient can be observed, as shown in Figure 10b. When the mass density increases, the damping coefficient increases at high frequency (ωh > s0). This implies that, at high frequency, if the mass density increases along the thickness direction, then the Lamb waves in the FGVM thin film should attenuate more quickly than those in a homogenous material. Conversely, if the Lamé parameters increase along the thickness direction, then the attenuation of Lamb waves in the FGVM thin film should be weakened. If only the relative viscosity coefficient increases along the thickness direction, then the attenuation of Lamb waves will become more serious, as shown in Figure 11. As the gradient coefficient increases, the damping coefficient increases and the attenuation tendency becomes obvious.

Conclusion
The power series method can be employed for solving the governing differential equations for Lamb wave propagation in FGVM thin films. The MMA method is proposed to solve the transcendental equations in the form of a series with complex coefficients. It is suggested that the meshing number should be selected as 6-8 in the first loop step and 4 in other loops. The numerical results obtained by these methods agree well with the exact analytical solution.
When Lamé parameters and mass density vary along the thickness direction identically, the influence of the gradient properties on the wave velocity is slight but that on the wave structure and the ellipticity of particles on the surface is obvious. This suggests that the ellipticity of particles on the surface should be selected to measure the gradient property if the bulk wave velocities are constants in the FGVM thin film. When Lamé parameters and mass density vary along the thickness direction independently, the variation of the phase velocity can be used for testing the gradient parameters. However, when the relative viscosity coefficient is a variable and both Lamé parameters and mass density are constants, the gradient property will not affect the phase velocity. The attenuation tendency becomes obvious with the increase of the gradient relative viscosity coefficient.

Conclusion
The power series method can be employed for solving the governing differential equations for Lamb wave propagation in FGVM thin films. The MMA method is proposed to solve the transcendental equations in the form of a series with complex coefficients. It is suggested that the meshing number should be selected as 6-8 in the first loop step and 4 in other loops. The numerical results obtained by these methods agree well with the exact analytical solution.
When Lamé parameters and mass density vary along the thickness direction identically, the influence of the gradient properties on the wave velocity is slight but that on the wave structure and the ellipticity of particles on the surface is obvious. This suggests that the ellipticity of particles on the surface should be selected to measure the gradient property if the bulk wave velocities are constants in the FGVM thin film. When Lamé parameters and mass density vary along the thickness direction independently, the variation of the phase velocity can be used for testing the gradient parameters. However, when the relative viscosity coefficient is a variable and both Lamé parameters and mass density are constants, the gradient property will not affect the phase velocity. The attenuation tendency becomes obvious with the increase of the gradient relative viscosity coefficient.
The method proposed herein and the results obtained should provide theoretical guidance for Figure 11. Damping coefficient of Lamb waves with the relative viscosity coefficient varying.

Conclusions
The power series method can be employed for solving the governing differential equations for Lamb wave propagation in FGVM thin films. The MMA method is proposed to solve the transcendental equations in the form of a series with complex coefficients. It is suggested that the meshing number should be selected as 6-8 in the first loop step and 4 in other loops. The numerical results obtained by these methods agree well with the exact analytical solution.
When Lamé parameters and mass density vary along the thickness direction identically, the influence of the gradient properties on the wave velocity is slight but that on the wave structure and the ellipticity of particles on the surface is obvious. This suggests that the ellipticity of particles on the surface should be selected to measure the gradient property if the bulk wave velocities are constants in the FGVM thin film. When Lamé parameters and mass density vary along the thickness direction independently, the variation of the phase velocity can be used for testing the gradient parameters.
However, when the relative viscosity coefficient is a variable and both Lamé parameters and mass density are constants, the gradient property will not affect the phase velocity. The attenuation tendency becomes obvious with the increase of the gradient relative viscosity coefficient.
The method proposed herein and the results obtained should provide theoretical guidance for ultrasonic nondestructive testing of heterogeneous viscoelastic materials and enable the safe evaluation of surface acoustic wave devices.