Application of Multi ‐ Parameter Perturbation Method to Functionally ‐ Graded, Thin, Circular Piezoelectric Plates

: In this study, a multi ‐ parameter perturbation method is used for the solution of a functionally ‐ graded, thin, circular piezoelectric plate. First, by assuming that elastic, piezoelectric, and dielectric coefficients of the functionally ‐ graded materials vary in the form of the same exponential function, the basic equation expressed in terms of two stress functions and one electrical potential function are established in cylindrical coordinate system. Three piezoelectric coefficients are selected as perturbation parameters, and the established equations are solved by the multi ‐ parameter perturbation method, thus obtaining up to first ‐ order perturbation solutions. The validity of the perturbation solution obtained is verified by numerical simulations, based on layer ‐ wise theory. The perturbation process indicates that adopting three piezoelectric coefficients as perturbation parameters follows the basic idea of perturbation theory—i.e., if the piezoelectricity may be regarded as a kind of introduced disturbance, the zero ‐ order solution of the disturbance system corresponds exactly to the solution of functionally ‐ graded plates without piezoelectricity. The result also indicates that the deformation magnitude of piezoelectric plates is smaller than that of plates without piezoelectricity, due to the well ‐ known piezoelectric stiffening effect.


Introduction
Piezoelectric materials have electromechanical coupling characteristics-i.e., they may generate mechanical deformation in an electric field and at the same time, electrical polarization under mechanical loads. This important characteristic makes them a good candidate for a variety of electromechanical devices-for example, sensors and actuators used extensively in electromechanical conversion [1][2][3][4]. Piezoelectric sensors are usually a laminated original made by ceramic slice. For this kind of laminated original, however, it is easy to cause stress concentration and promote the growth of interfacial microcracks, which further limits the application and development of the piezoelectric original. In order to solve this problem, functionally-graded piezoelectric materials (FGPMs), whose material properties change continuously in one direction, have been developed [5][6][7]. Since there is no obvious interface in this kind of material, the damage caused by the stress concentration at the interface can be avoided. This undoubtedly promotes the application of piezoelectric materials in electromechanical conversion.
With the increasing application of functionally-graded piezoelectric materials, precise characterization of their mechanical properties is urgently needed. A great deal of research has been done on the mechanical properties of functionally-graded piezoelectric materials and structures: for example, FGPM cantilever beams [8][9][10][11][12], FGPM plates [13][14][15][16][17], and FGPM shells [18][19][20][21]. At the same time, the generation of new problems also puts forward greater requirements for the corresponding solving methods. Among all kinds of solving methods for the problem (analytical, numerical, even experimental, for example), the analytical approach is relatively challenging for the complexity of the governing equations obtained. It is known that due to electromechanical coupling, the constitutive equation becomes more complicated, compared with the purely mechanical problem in classical elasticity. The governing equation thus obtained and the variation of material properties along certain directions, due to the incorporation of functionally-graded materials, make analytical solving more and more difficult. It has been found that even for a plate consisting of classical materials, the analytical solution is hard to obtain, and in most cases, numerical simulations have to be resorted to (for example, see recent studies [22,23]). Therefore, seeking an effective analytical technique, especially since the structural property of the analytical solution may clearly show the electromechanical coupling effect in the final results, seems to be valuable and urgent. The purpose of this paper is based on this key point.
The perturbation method is one of the standard analytical methods used for the solution of nonlinear problems in applied mechanics and physics. This method consists of developing the solution of an initial or boundary value problem in an asymptotic series of a parameter, which appears explicitly either in the problem or is introduced artificially. More specifically, during the perturbation, the unknown functions are expanded in the form of ascending powers with respect to a certain small parameter. Substituting these expansions into the governing equations and corresponding boundary conditions will yield a series of equations used for determining the approximate solution of all levels by equating the same order of the perturbation parameter.
Although perturbation solutions may not rely on any small parameter, as suggested by Van Gorder [24], the perturbation parameter plays an important role in perturbation, since the right choice may permit us to obtain asymptotic solutions with better convergence. The earlier classical works may be found in Vincent's and Chien's studies. By using the external load as a perturbation parameter, Vincent [25] first obtained the perturbation solution of the Föppl-von Kármán equations. Given that the perturbation parameter either appears explicitly or is introduced artificially into the problem, Chien [26] obtained another perturbation solution using the central deflection as a perturbation parameter. Compared with experimental results, Chien's solution is accurate, and it has been cited as a classical work in subsequent studies for a long period of time. In addition to the load and the central deflection, there are several other choices for perturbation parameters-for example, a generalized displacement [27], a linear function of Poisson's ratio [28], and an average angular deflection [29]. Chen and Kuang [30] discussed the differences between the possible perturbation parameters.
In the face of the confusion of parameter selection, we usually have two entirely different ways of solving it. One is the non-parametric perturbation method proposed by Chen [31,32], in which the physical meaning of the small parameter is not given in advance, thus avoiding the empirical factors in the process of parameter selection. In other words, the selection of parameters is not constrained and is thus free, so it is also referred to as the free parameter perturbation method. Another method is the so-called "multi-parameter perturbation" method, in which two or more parameters are selected as perturbation parameters. The earlier work in this method may be found from Nowinski and Ismail [33], in which Föppl-von Kármán equations of elastic anisotropic plates were solved by selecting the anisotropy of materials and the load as two perturbation parameters. The first application of the biparametric perturbation method to beam problems was performed by Chien [34], in which the load and the height difference of end supports were selected as two perturbation parameters, and thus the classical Euler-Bernoulli equation was solved. Later, He and Chen [35] derived a biparametric perturbation solution for the same problem by simplifying the governing equation. He et al. [36] further studied the application of the biparametric perturbation method to beam problems with the height difference of end supports under various boundary conditions. More recently, progress has been made in the application of this multi-parameter perturbation to bimodular plates, mainly including the combined loads problem [37], the parameter selection problem in a single load [38], and the parameter selection problem in a functionally-graded material (FGM) plate [39]. More recently, Fallah et al. [40] reviewed the perturbation method in mechanical, thermal, and thermo-mechanical loadings of cylindrical bendings of FGM plates, both clamped and simply supported, in which one-and two-parameter perturbations were used. It is worth pointing out that Lian et al. [41] used a three-parameter perturbation method to solve a FGPM cantilever beam problem for the first time, in which three piezoelectric coefficients, 31 D , 33 D , and 15 D , were selected as perturbation parameters. This may be the first application of a multi-parameter perturbation in FGPM beams. Adopting piezoelectric coefficients as perturbation parameters follows the basic idea of perturbation theory, which has been reported in our previous study [42]. This basic idea indicates that so-called perturbation is essentially a kind of disturbance for an undisturbed system: if the solution of the undisturbed system may be obtained, the corresponding solution of that disturbed system may also be obtained via the perturbation technique. Given that there has been no similar application in FGPM plates, this study seems to be necessary.
In this study, we will use the multi-parameter perturbation method to solve the bending problem of a functionally-graded, thin, circular, piezoelectric plate. This paper is organized as follows. First, in Section 2, the basic equation expressed in terms of stress functions and potential function is established, in which elastic, piezoelectric, and dielectric coefficients of the functionally-graded materials vary in the form of same exponential function. By selecting three piezoelectric coefficients as perturbation parameters, in Section 3, we obtain up to first-order perturbation solution under boundary conditions. Numerical simulation and comparison with perturbation solution are given in Section 4. The effect of the gradient index on the solution and deflections of FGPM and FGM plates are discussed in Section 5. Three important conclusions and subsequent works are summarized in Section 6.

Mechanical Model and Basic Equations
In this paper, we will study the bending problem of a functionally-graded, thin, circular, piezoelectric plate with radius a and thickness h, as shown in Figure 1, in which the thin circular plate is subjected to uniformly-distributed external load q on its upper surface, and in which the boundary conditions of its peripheries are considered to be fully fixed. Obviously, an axisymmetric deformation will generate in this plate, and it is convenient to introduce the cylindrical coordinate system to describe this problem; thus, an r z    coordinate system is established, in which the coordinate origin o is set at the geometric center of the plate, r is the radial coordinate,  is the circumferential coordinate (which is not shown due to axisymmetric feature), and z stands for the vertical direction, as shown in Figure 1.

Nondimensionalization and Perturbation Expansions
We introduce the following dimensionless quantities: where Equation (15a) is used for the nondimensionalization of governing equations, and Equation (15b) is only for later computation. Substituting the above dimensionless quantities into the three governing equations, (9), (10), and (11), and for the convenience of writing and reading, also replacing these dimensionless quantities 33 , , , , , , , , , , , , , , with the original forms, respectively, 33 , , , , , , , , , , , , , ,  ( 2 ) 4  4  2  2  3  3  23  3  3  13  12  31  15  31  4  2  2  2  2   3  3  3  3  15  33  31  2 3 and  2  3  3  2  3  2  2  2  13  13  13  12  2  3  2  2  3   2  2  31 (1 ) In the above three equations, , , and are dimensionless constants. If there is no special explanation below, the next symbols in this paper refer to dimensionless quantities. We select the zero point values of three piezoelectric coefficients, , , and , as our perturbation parameters, and expand the unknown functions, , , and , with respect to the parameters: , , We also let where , , , , and , are the j-th functions with respect to the i-th order of perturbation parameters; and , , , , , ℎ , ℎ , and ℎ are the component functions. It should be noted here that the expansion pattern of , , , , and , -i.e., Equation (20)-is consistent with the pattern of a similar beam problem [42]. In addition, we found that the odd power term of  disappears in the later derivation, and only even power terms remain. For the stress function ( , )    , 0  term will disappear in the later stress derivations (easily seen from Equation (8)); thus, there is no need to list this term in the expansion of ( , )    .
The boundary conditions used for the zero-order perturbation solution are, for the upper surface: for the lower surface: and for the fixed peripheral: 0, 0, 0, at and .
The integration constants are determined as Note that (0) 5 C is the constant term of electrical potential function  , and has no effect on the electrical field (this may be easily seen from Equation (6)); therefore there is no need to determine it. This conclusion is also applicable to the next first-order perturbation. In addition, we also note that unlike Equations (31)-(32), Equation (33) is not expressed in terms of two stress functions. It may be more convenient to use this boundary condition directly expressed in terms of displacement. For example, (0)

5
A and (0) 6 A are still not determined, and will be given in Appendix A, in which Equation (33) is used.

Numerical Simulation and Comparison with Perturbation Solution
In order to further verify the correctness of the multi-parameter perturbation solution, a FGPM thin circular plate subjected to uniformly-distributed loads was calculated via a numerical analysis program, thus allowing corresponding comparison with the analytical expression obtained above.

Numerical Simulation
ABAQUS software is one of the large-scale finite element software programs at present that can analyze complex mechanical problems in engineering, including the problem of piezoelectric materials. However, the software itself does not include functionally-graded properties of materials varying in a form of continuous function. For this purpose, we resorted to layer-wise theory to simulate the problem studied.
Another problem, which should be noted here, is different types of constitutive models of the piezoelectrical materials. In ABAQUS, an e-type constitutive equation for piezoelectrical materials is adopted, such that note that in the theoretical derivation presented above and the properties of materials given below, the d-type constitutive equation of piezoelectric materials was adopted. Therefore, the transformation of these material parameters in different types of the constitutive model was needed. Considering the length limit of the paper, however, we will not repeat it here. More details may be found in our previous study [45]. The other detailed steps for modeling and simulation are as follows: (i) Establishment of entity structure The solid model of an FGPM thin circular plate is established, in which the thickness of the plate 20mm h  and the radius of the plate 300mm a  .

(ii) Layering and determination of materials properties in each layer
According to conventional practice, we still adopted layer-wise theory to simulate the functionally-graded properties varying in the thickness direction. Without losing computational accuracy, the plate was divided into a moderate number of layers; the physical parameters of the material in each layer was regarded as the same, thus indirectly realizing the continuous variation of material properties throughout the thickness direction. The change in material properties in each layer will follow the form that has been adopted in our theoretical derivation, where the functionally gradient index is set to 2   . The material constant at 0   is shown in Table 1, which may be input into ABAQUS directly.  In our real simulation, a plate with thickness 20 mm is divided into 19 layers along the thickness direction, in which the thickness of the middle layer is 2 mm, the upper and lower layer numbers are both 9 mm, and the thickness of each layer is 1mm, as shown in Figure 2. Therefore, via the layered model, the physical parameters in each layer may be input into the program layer by layer.

(iii) Establishment of boundary conditions
The boundary condition of the peripheral of the plate was still considered as fully fixed, as shown in Figure 1. For this purpose, we needed to exert constraint at the periphery of the plate, including displacement and rotation.

(iv) Mesh division
In this simulation, an eight-node, linear piezoelectric brick C3D8E was adopted. Considering that the assumed material property varies more uniformly, mesh of the same density was adopted. The global approximate size was set to be 0.03.
(v) Step module and adding loads An analysis step called ʺstatic general" was established to apply the load. In this simulation, only a uniformly distributed load form was considered. For this purpose, a uniformly distributed load along the positive direction of z axis was applied on the upper surface of the plate.

(vi) Operation and results output
After a job was established, the job was submitted for calculation and output of the results.

Comparison with Perturbation Solution
In this section, we will use the results obtained in Section 3 to compute the theoretical solution for some special quantities in question, in which the gradient index 2   and other given data is the same as numerical simulation (for example, the physical properties of PbZrTiO3-4 (Generally abbreviated as PZT-4) in a dimensionless form are given in Table 2, according to Equation (15a)), thus comparing the theoretical solution with the numerical example in a dimensionless form, as shown in Figure 3. Note that all stress components in Figure 3 are given in dimensionless values (please refer to Equation (15b)).  It is easy to see that the theoretical solution and numerical result increase nonlinearly with the distance from the plate edge; in particular, the two deflection values achieve their maximum at the plate's center. The two curves agree well with each other, and the theoretical curve is slightly higher than the numerical one. Figure 3b shows the variation of the radial stress r  at the plate center 0   , with thickness direction coordinate  . The two curves basically agree, and there are positive and negative radial stresses in them. The location at which the radial stress is zero is exactly the location of the neutral layer. Obviously, the neutral layer is not at 0   . The greater the distance from the neutral layer, the greater the radial stress. This is similar to the distribution of radial stress in general bending problems.  ). The two curves have the same variation tendency; however, there is still some difference between them. It is easy to see that in the area near the upper surface, or at the load action surface, the numerical results are quite different from the theoretical values, which is mainly due to the large discreteness of the finite element calculation results near the load action surface. Fortunately, z-direction stress in thin plates is a relatively secondary stress, like the extrusion stress along the vertical direction in the beam. Compared with the radial stress and circumferential stress, the value of z-direction stress is negligibly small; this may be seen by observing the values of the coordinate axis from Figure 3b-d. In addition, Table 3 gives the maximum values of vertical deflections, radial and circumferential stresses, and the relative errors between theoretical solutions and numerical results. It can be seen from Table 3 that the error of deflection is very small, and the error of the maximum circumferential stress is relatively large, which is also related to the large discreteness of the numerical results of the finite element method on the upper surface on which the load acts.

Results and Discussions
After the validity of the multi-parameter perturbation solution is verified by numerical simulation, the perturbation solution may be used to discuss some interesting topics, including the effect of the functional gradient index on the solution, as well as the deflection difference between FGPM plates and FGM plates.

Effect of Gradient Index on Solution
To analyze the effect of the gradient index on the solution, we select some values of gradient index, i.e., 2, 1, 0.01, 2, and 5, in which 0.01   may approximately stand for the case that the gradient index is zero, since  cannot be set to zero in our numerical computation. For this purpose, the variations of mechanical stresses and displacements, as well as electrical potentials, electrical field intensities, and electrical displacements along the direction of radius or thickness, are plotted in different gradient indices, as shown in surfaces, and achieves the greatest value at the neutral layer, which is consistent with classical problems. Figure 4d shows the variation of the deflection z u at the neutral layer along the direction of radius. It is easy to see that under various gradient indices, the maximum deflection is also different. Obviously, due to the difference introduced by the physical parameters of functionally-graded material, the bending stiffness of the plate is different for various gradient indices.  field intensity at the plate's edge, (f) electrical potential at the lower surface, (g) radial electrical displacement at the plate's edge, and (h) z-direction electrical displacement at middle layer.
From Figure 4e, it can be seen that the distribution of r E at the plate's edge under various gradient indices, along the direction of thickness, is basically similar, only with a slight change along the thickness direction. Figure 4f shows the distribution of  at the lower surface along the direction of the radius. Due to the connection to a ground, the electrical potential is zero at the plate's edge ( 15   ) and increases with the distance from the plate edge to the plate center, finally achieving the largest value at the plate center ( 0   ), with the exception of the case where 2 in which the maximum value is obtained at about 7.5   . Figure 4g shows that the radial electrical displacement at the plate edge varies with the vertical coordinate. r D changes monotonously along the direction of thickness. When the gradient index is positive or negative, the tendency for change is contrary. Figure 4h shows the variation of the z-direction electrical displacement at the middle layer with the radial coordinate. It can be seen that under various functional gradient indices, z D at the plate center is zero, and increases gradually from the center to the edge of the plate, achieving its maximum value at the edge. We note that although r D and z D are attributable to electric displacement components, and may be easily found, as far as their magnitude is concerned. This fact may be used for the simplification in some cases.

Deflection Comparison Between FGPM and FGM Plates
From the multi-parameter perturbation solution for the three piezoelectric coefficients 31 D , 33 D , and 15 D , we may easily obtain the solution of purely functionally-graded thin plates without the piezoelectric effect, by only letting the three parameters be zero-or, alternatively, by only taking the zero-order perturbation solution in Equation (19). Thus, we have Due to the fact that the perturbation solution is, in this form, an expansion of the power series with respect to the perturbation parameter, the unperturbed solution (zero-order perturbation solution) may be approximately regarded as a linear part to the solution of the perturbation system (first-order up to the higher-order perturbation solution). Here we focus our discussion on deflection comparisons between FGPM plates and FGM plates; thus, locations of the neutral layer and central deflections of FGPM plates and FGM plates, under various functionally gradient indices, are computed and listed in Table 4. It is easy to see that under various gradient indices, the locations for the unknown neutral layer are also different: the greater the absolute value of the gradient index, the greater the distance from the geometric center layer of the plate. In particular, when the absolute values of the gradient indices are equal ( 2    and 2   , for example), the position of the neutral layer is symmetrical about the geometric middle layer. This conclusion has been obtained from our above discussion on Figure 4a-c. Due to the fact that, under the influence of functionally-graded materials, the neutral layer no longer coincides with the geometric middle plane, the position of the neutral layer is generally close to the layer with a large stiffness coefficient. In addition, it is easy to see that under various gradient indices, there is no more difference in central deflection values for FGM plates, while for FGPM plates, the change of the central deflection is relatively large.
For further analysis, we may let the deflection at the neutral layer, z u , be the following form, according to our obtained solution: where a, b, and c are constants, and their values under various gradient indices are listed in Table 5. It is obvious that the constant a corresponds exactly to the central deflection values in Table 4, which is easily seen from Equation (54).  In addition, we may plot the deflection curves of FGPM plates and FGM plates under various functional gradient indices, as shown in Figure  FGM(1) are also pretty close. For the FGPM plates, this is not the case. This indicates that the change of functional gradient index may influence the deflection of FGPM plates to some extent, while for FGM plates, there seem to be no more differences when the gradient index varies. It is worth noting that under the same functional gradient index, the deflection value of FGPM plates is smaller than that of corresponding FGM plates, which shows that the intervention of the piezoelectric effect may reduce the deformation of the plate. This phenomenon may be further explained from the point of view of energy conservation and transformation. For FGPM plates, a portion of the work done by applied external loads is transformed into electrical energy, due to the piezoelectric effect of materials; therefore, that portion transformed into elastic strain energy decreases correspondingly, and the deformation magnitude will be smaller than that of FGM plates, in which there is no so-called electrical energy transformation.

Conclusions
In this study, we use a multi-parameter perturbation method to solve a functionally-graded, thin, circular piezoelectric plate under the action of uniformly-distributed loads, in which three piezoelectric coefficients, 31 D , 33 D , and 15 D , are selected as perturbation parameters. The validity of the multi-parameter perturbation solution obtained was verified by numerical simulations based on layer-wise theory. The following three conclusions can be drawn: (i) Adopting the three piezoelectric coefficients as perturbation parameters follows the basic idea of perturbation theory-i.e., if pure FGM plates without piezoelectric effects are taken as the undisturbed system, and the piezoelectric effect may be introduced as a kind of disturbance, then FGPM plates may be regarded as a disturbed system; thus, the solution for pure FGM plates may be easily obtained as a zero-order solution of the disturbed system. (ii) In our perturbation, two stress functions and one electrical potential function were selected as basic functions. It was found that in the zero-order perturbation solution, only stress functions were not zero; thus, the zero-order solution actually corresponded to the elastic solution concerning elastic stress, elastic strain, and elastic displacement, this conclusion is consistent with the former conclusion: while in the first-order perturbation solution, only the electrical potential function was not zero; thus, the first-order solution actually corresponded to the electrical solution, concerning electrical potential, electrical field intensity, and electrical displacement. (iii) The deformation magnitude of FGPM plates is generally smaller than that of corresponding FGM plates, showing the well-known piezoelectric stiffening effect. From the point of view of energy conservation and transformation, a portion of the work done by applied external loads is transformed into electrical energy, due to the piezoelectric effect, thus decreasing elastic strain energy stored in FGPM plates and resulting in the corresponding decrease in deformation magnitude.
It is undeniable that the calculation process of the multi-parameter perturbation method is somewhat cumbersome; this is the reason why we only calculate up to the first-order perturbation solution. If a higher accuracy of the solution is required, we need to resort to a higher order perturbation solution. However, the multi-parameter perturbation method proposed in this study may provide a reference for solving multi-physical field problems like electric, magnetic, or thermal fields, in addition to the traditional mechanical field. The theoretical result presented in this study is helpful to precisely analyze the mechanical properties of functionally-graded piezoelectric materials and structures, as well as to design sensors and actuators used extensively in electromechanical conversion. The related work is in progress.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results. 16 (1 ) 16 (1 ) 12 4 1 the zero-order electrical displacements and electrical strength solutions 15 31 31 33 as well as the radial displacement solution (1 ) (1 ) Note that (0) 5 A and (0) 6 A , as well as z-direction displacement z u , are still not determined.
Next, we will derive the numerical value solution of (0 (A5) Lastly, we obtain which may be used for the solution of unknown neutral layer n  . Table A1 gives values of n  in different gradient indexes. A in different gradient indexes, as shown in Table A2.
Note that there is no need to determine (1) 5 C , since it is the constant term of electrical potential function and has no effect on the electrical field.
Similarly, for D15, substituting Equation (20) into Equation (45) Note that there is also no need to determine (1) 17 C .