One-Dimensional Theoretical Solution and Two-Dimensional Numerical Simulation for Functionally-Graded Piezoelectric Cantilever Beams with Different Properties in Tension and Compression

The existing studies indicate polymers will present obviously different properties in tension and compression (bimodular effect) which is generally ignored because of the complexity of the analysis. In this study, a functionally graded piezoelectric cantilever beam with bimodular effect was investigated via analytical and numerical methods, respectively, in which a one-dimensional theoretical solution was derived by neglecting some unimportant factors and a two-dimensional numerical simulation was performed based on the model of tension-compression subarea. A full comparison was made to show the rationality of one-dimensional theoretical solution and two-dimensional numerical simulation. The result indicates that the layered model of tension-compression subarea also makes it possible to use numerical technique to simulate the problem of functionally graded piezoelectric cantilever beam with bimodular effect. Besides, the modulus of elasticity E* and the bending stiffness D* proposed in the one-dimensional problem may succinctly describe the piezoelectric effect on the classical mechanical problem without electromechanical coupling, which shows the advantages of one-dimensional solution in engineering applications, especially in the analysis and design of energy harvesting/sensing/actuating devices made of piezoelectric polymers whose bimodular effect is relatively obvious.


Introduction
Piezoelectric materials have an electromechanical coupling characteristic, which makes them a good candidate for a variety of electromechanical devices, for example, sensors and actuators used extensively in electromechanical conversion. Piezoelectric sensors are usually a laminated original made by ceramic slice, so it is easy to cause stress concentration and promote the growth of interfacial microcracks. In order to overcome this difficulty, scholars developed functionally graded piezoelectric materials (FGPM) whose properties of materials change continuously along certain direction. There is no obvious interface in FGPM, thus the damage due to the stress concentration at the interface is effectively avoided. Studies on FGPM have attracted the attention of scholars from all over the world. In recent years, with the development, universality and miniaturization of electronic devices, new piezoelectric materials continue to emerge, among which piezoelectric polymers play an increasingly important role [1][2][3][4]. Piezoelectric polymers are attractive for wearable due to their flexibility and with different moduli in tension and compression. Although the effectiveness of the analytical work was verified by comparing with the existing theoretical work presented by Zhong and Yu [13], it is also a pity that the comparison was based on the degraded analytic expressions from He et al. [27] since the work of Zhong and Yu [13] did not consider the bimodular effect. Now that existing analytical works are not fully qualified for this comparative work, we have to resort to the numerical technique based on existing software such as ABAQUS. For a traditional piezoelectric problem, the use of ABAQUS appears to be only a step-by-step process. However, how about after introducing the bimodular effect and functionally graded properties of the materials? In addition, similar to the bending problem of classical beams, a so-called one-dimensional solution and two-dimensional solution under different application conditions always exists. Now that the two-dimensional theoretical solution for a bimodular FGPM beam can be obtained [27], what is the form of one-dimensional solution? For the above two reasons, we think it is necessary to obtain the one-dimensional theoretical solution and also perform the two-dimensional numerical simulation for the problem of bimodular FGPM cantilever beam. Therefore, this study may serve as a new supplement to the existing works, not only from the theoretical aspect but also from the point of view of numerical simulation.
In this study, we will derive one-dimensional theoretical solution by neglecting some unimportant factors and perform two-dimensional numerical simulation based on the model of tension-compression subarea. The whole paper is organized as follows. In Section 2, the solving problem will be described, including the definition of different properties in tension and compression and the constitutive equation of FGPM in a two-dimensional case. The one-dimensional theoretical solution will be derived in Section 3 and the two-dimensional numerical simulation will be performed in Section 4. Next, in Section 5, we will make extensive comparisons with existing studies which not only include previous studies of our themselves [27] but also the work from other authors, to show the validity of our work, and also study the evolution from classical beams to bimodular FGPM beams as well as discuss the deformation of flexible piezoelectric structures. According to the results allude to above, some main conclusions will be drawn in Section 6.

The Problem Description
As indicated above, the characteristics of material considered in this study include the bimodular effect, the functionally graded property and the piezoelectric characteristic. In practical applications, there is a large number of materials containing three properties above, for example, the functionally graded material made of piezoelectric ceramics and steel, at the same time, the bimodular effect is considered due to the existence of ceramics, which present relatively obvious different elastic properties in tension and compression (see [23]). Among these applications, they usually exist in a certain structural form, for example, the form of cantilever beam; thereby it is necessary to study a bimodular functionally graded piezoelectric cantilever beam.
An orthotropic functionally graded piezoelectric cantilever beam with bimodular effect is considered here, as shown in Figure 1, in which the right end of the beam is fully fixed and the left end free; h × b stands for the rectangular section dimension of the beam and l is the length of the beam (h l). Without losing generality, the loads may be considered as single form load, for example, the upper layer of the beam is subjected to uniformly distributed load q, or the left end of the beam is subjected to a shear force P or a bending moment M, or the combined loads from above load forms, as shown in Figure 1. It is obvious that the loads acting in plane coordinate system xoz may cause downward bending of the beam, thereby generating the so-called tensile part and compressive part, bounded by the neutral layer. Thus, we establish a rectangular coordinate system in which z = 0 is exactly at the neutral layer, as shown in Figure 1. The upper and lower edge layers are z = −h 2 and z = h 1 , respectively, in which h 1 is the tensile height and h 2 the compressive height according to previous studies [28,29]. Note that due to the introduction of functionally graded property, physical parameters of materials of the beam are also functions of coordinates. Generally, it is assumed that physical parameters vary only along certain direction, for example, the thickness direction. In this study, material parameters are assumed to vary with z, in the light of the following rules where F + (z) = e α 1 z/h , F − (z) = e α 2 z/h are gradient functions in tension and compression, respectively; superscript "+" denotes a tensile quantity and "−" compressive quantity; ij are values of corresponding material parameters at the neutral layer (z = 0), respectively. Note that a set of very small electrodes are adhered discontinuously to the upper and lower surfaces of the beam, and at the same time the beam is poled along the direction of z. In two-dimensional problem, let the stress components and strain components be σ +/− zx , respectively; let the electrical displacement components and the electrical field components be D +/− respectively, thus the physical equations give and where superscript "+/−" still denotes a tensile (compressive) quantity, similar to Equation (1). It may be inferred that the constitutive relation in two-dimensional case may be moderately simplified in one-dimensional problem, according to our previous study concerning one-dimensional and two-dimensional problems of bimodular FGM beams [30].
direction. In this study, material parameters are assumed to vary with z, in the light of the following rules and where superscript "+/−" still denotes a tensile (compressive) quantity, similar to Equation (1). It may be inferred that the constitutive relation in two-dimensional case may be moderately simplified in one-dimensional problem, according to our previous study concerning one-dimensional and twodimensional problems of bimodular FGM beams [30].

One-Dimensional Theoretical Solution
In one-dimensional problem, we may consider the simplest case, i.e., the pure bending problem shown in Figure 2, in which the left end of the cantilever beam is subjected to a bending moment M. Obviously, the beam will generate downward bending under the action of the bending moment in plane coordinate system xoz, thus forming a so-called tensile zone and compressive zone. We still establish the neutral layer at 0 z = , and the tensile modulus of elasticity is denoted by

One-Dimensional Theoretical Solution
In one-dimensional problem, we may consider the simplest case, i.e., the pure bending problem shown in Figure 2, in which the left end of the cantilever beam is subjected to a bending moment M. Obviously, the beam will generate downward bending under the action of the bending moment in plane coordinate system xoz, thus forming a so-called tensile zone and compressive zone. We still establish the neutral layer at z = 0, and the tensile modulus of elasticity is denoted by E + (z) and the compressive one by E − (z); similarly, the tensile and compressive heights are h 1 and h 2 , respectively, as shown in Figure 2.

Mechanical Stress and Deflection
Note that in such a one-dimensional pure-bending problem, there is no stresses In existing studies for two-dimensional problem [13,14,27] Substituting Equation (6) into Equation (4), we obtain where * E is defined as the modulus of elasticity in one-dimensional problem, i.e., , which is exactly the reciprocal relationship between flexibility coefficient and stiffness coefficient. From the viewpoint of regression satisfaction, this fact verifies indirectly the correctness of Equations. (4-5) in one-dimensional case. After substituting the functionally graded form, i.e., Equation (1), into Equation (7), we have

Mechanical Stress and Deflection
Note that in such a one-dimensional pure-bending problem, there is no stresses σ +/−  (2) and (3), may be simplified as In existing studies for two-dimensional problem [13,14,27], D x >> D z may be found, thus we may assume D z ≈ 0 in the one-dimensional problem. From Equation (5), we have Substituting Equation (6) into Equation (4), we obtain where E * is defined as the modulus of elasticity in one-dimensional problem, i.e., We note that if the above equation is rewritten as the form , it may clearly explain the piezoelectric effect on the modulus of elasticity in classical problem. Specially, when

33
(for example, for the piezoelectric materials PZT-4, it is the case), we may have E * = 1/s +/− 11 , which is exactly the reciprocal relationship between flexibility coefficient and stiffness coefficient. From the viewpoint of regression satisfaction, this fact verifies indirectly the correctness of Equations (4) and (5) in one-dimensional case. After substituting the functionally graded form, i.e., Equation (1), into Equation (7), we have where e α i z/h are gradient functions, introduced earlier in Equation (1), α i (i = 1, 2) corresponds to "+/−"; and s 0 11 ,λ 0 33 and d 0 31 represent the values of material parameters s 11 ,λ 33 and d 31 at the neutral layer (z = 0), respectively.
If we let the curvature radius of the beam in bending be ρ, the strain at any point may be expressed as, according to plane section assumption in a one-dimensional case, Substituting Equation (10) into Equation (9), we obtain Thus, we obtain the bending stress in tensile and compressive zone for the one-dimensional pure-bending problem, i.e., for 0 ≤ z ≤ h 1 , and for −h 2 ≤ z ≤ 0, Note that ρ, h 1 and h 2 , as well as the deflection of the beam are still not determined yet. Next, we will use the conditions of internal forces on the section to determine them. Let the normal internal force acting on any section be N, thus N = 0 will give Substituting Equations (12) and (13) into Equation (14), we have where k 0 b/ρ is a constant and may be deleted in the above equation. If we let from Equation (15), we may obtain which is used for solving the tensile and compression section height, i.e., h 1 and h 2 .
The bending moment acting on any section is M(x) = M, this will give We have, after substituting Equations (12) and (13) into Equation (18) thus from Equation (19), we have a familiar form concerning deformation as follows where D * is defined as the bending stiffness of a bimodular FGPM beam, . Substituting Equation (21) into Equations (12) and (13), we have, for 0 ≤ z ≤ h 1 , If we let the deflection be w, Euler-Bernoulli equation in small-deflection case will give Integrating the above equation with respect to x will yield (M is a constant) where c and d are two integrating constants and may be determined by the following boundary conditions: Thus, Substituting it into Equation (25), we finally obtain For the convenience of the next comparison, when the cantilever beam is subjected a uniformly-distributed load on its upper surface, M(x) = qx 2 /2, the two integrating constants may be again determined as Thus, we have w(x) = q Polymers 2019, 11, 1728 8 of 24

Electrical Displacement
Next, we will derive the electrical displacement components D +/− x in one-dimensional case. Note that the electrical displacement is generated not only from the electrical voltage application in piezoelectric elements, as indicated in actuator model, but also from the mechanical load, as indicated in sensor model which agrees with our study model, according to Figure 2 (under the action of bending moment).
From Equations (1), (6), (22) and (23), we have where l 0 = d 0 31 /λ 0 33 . Let the potential function be Φ +/− , the relation of electrical field and potential in two-dimensional problem will give is an unknown function concerning only x. However, in the one-dimensional problem, we think the variation of Φ +/− with x is only embodied in the bending moment which has been included in E +/− z (see Equation (31)). This fact may be further demonstrated, from the side, based on the previous work [31] concerning purely piezoelectric materials without bimodular functionally graded properties, in which the potential function was determined as, in a two-dimensional case [31], where I is the moment of inertia of cross section, P is the concentrated force acting on the left end of the beam thus Px stands for the bending moment. Obviously, Substituting Equation (34) into the first expression of Equation (32) and also noting M(x) = qx 2 /2 if uniformly-distributed load is still considered here, E +/− x may be determined as, for 0 ≤ z ≤ h 1 , and for −h 2 ≤ z ≤ 0, From the first expression of Equation (5), Finally, we obtain the bending stress σ +/− x , the vertical deflection w and electrical displacement D +/− x in the one-dimensional problem.

Two-Dimensional Numerical Simulation
In this section we will use the software ABAQUS to simulate a bimodular FGPM cantilever beam subjected to a uniformly-distributed load q on its upper surface, as shown in Figure 3, in which other physical quantities and the establishment of coordinate system are the same as those in Sections 2 and 3. At the same time, it is assumed that 1 corresponds to the x direction, 2 to the y direction and 3 to the z direction.
From the first expression of Equation (5) and for Finally, we obtain the bending stress / x σ + − , the vertical deflection w and electrical displacement / x D + − in the one-dimensional problem.

Two-Dimensional Numerical Simulation
In this section we will use the software ABAQUS to simulate a bimodular FGPM cantilever beam subjected to a uniformly-distributed load q on its upper surface, as shown in Figure 3, in which other physical quantities and the establishment of coordinate system are the same as those in Sections 2 and 3. At the same time, it is assumed that 1 corresponds to the x direction, 2 to the y direction and 3 to the z direction.

Constitutive Equation of Piezoelectrical Materials
In ABAQUS, e-form constitutive equation of piezoelectrical materials is adopted, such that The z axis is set to be the polarization direction, corresponding to the 3-direction of ABAQUS. Since the elastic properties of ABAQUS are represented by the stiffness coefficient matrix ijkl D of corresponding material, which is the inverse matrix of the flexibility matrix, it is necessary to transform the flexibility coefficient matrix of the above material into the stiffness coefficient matrix, i.e.,

Constitutive Equation of Piezoelectrical Materials
In ABAQUS, e-form constitutive equation of piezoelectrical materials is adopted, such that where σ ij is the stress component; ε ij is the strain component; q ij is the electrical displacement is the piezoelectrical stress constants matrix; e ϕ mij is the dielectric constant matrix; E m and E j are electrical field strength.
The z axis is set to be the polarization direction, corresponding to the 3-direction of ABAQUS. Since the elastic properties of ABAQUS are represented by the stiffness coefficient matrix D ijkl of corresponding material, which is the inverse matrix of the flexibility matrix, it is necessary to transform the flexibility coefficient matrix of the above material into the stiffness coefficient matrix, i.e., Note that since we define the variation form of functionally graded materials in advance, the flexibility coefficient s ij = s 0 ij e α i z/h , where i = 1, 2 due to the bimodular effect, thus, the stiffness coefficient will change as c ij = c 0 ij e −α i z/h , otherwise s ij and c ij cannot [E] is an unit matrix. At the same time, piezoelectrical strain constants where c E is short circuit elastic stiffness constant matrix, thus piezoelectrical stress constants matrix is The constitutive equation of piezoelectric materials is expressed as follows, in the form of matrix, In ABAQUS, the double-subscript second-order tensor mark, 11, 22, 33, 13, 23 and 12, correspond to the vector components, 1, 2, 3, 4, 5 and 6, respectively. Thus, the above two equations are expressed as, in ABAQUS and where, D ijkl is the modulus of elasticity, D ij is the dielectric coefficient and q i is the electrical displacement component. Through the comparison of the above two sets of equations, it is easy to see the corresponding relationship of constants, which can be used to input values of the constants. For example, we should input c 11 at the location of D 1111 in ABAQUS; input e 31 at the location of e 311 and input λ 11 at the location of D 11 .

Modeling and Simulation
ABAQUS software is one of the large-scale finite element software at present which can analyze complex engineering mechanics problems including the problem of piezoelectric materials. However, the software itself does not involve different properties in tension and compression and functionally graded properties of materials varying in a form of continuous function. For this purpose, we should resort to subareas model in tension-compression and layer-wise theory to simulate the problem studied, which inevitably complicates the analysis process further. The detailed steps for modeling and simulation are as follows.
(i) Establishment of entity structure The solid model of a FGPM cantilever beam is established, in which the length of the beam l is set to be 1 m, the section height h to be 0.2 m and the section width b to be 0.08 m.
(ii) Determination of neutral layer and tension-compression subarea Note that comparing to common piezoelectric cantilever beam, the subarea of tension and compression of beam under external loads is a newly introduced feature, we must determine the unknown neutral layer first, thus realizing the so-called tension-compression subarea. For this purpose, a set of functionally graded indexes should be chosen to determine the tensile height and the compressive one of the beams, according to the one-dimensional theoretical solution. For example, consider the case α 1 = −2 and α 2 = −3. Substituting the given values into Equations (16) and (17) and also noting the section total height h = 0.2 m, we may have h 1 = 0.06 m and h 2 = 0.14 m. It is easy to see that due to the tensile modulus E + (z) = E 0 e −α 1 z/h is totally greater than the compressive modulus E − (z) = E 0 e −α 2 z/h , the neutral layer will locate below the geometrical middle surface, i.e., h 1 < h 2 , see Figure 4, in which E 0 is the modulus of the neutral layer.

Modeling and Simulation
ABAQUS software is one of the large-scale finite element software at present which can analyze complex engineering mechanics problems including the problem of piezoelectric materials. However, the software itself does not involve different properties in tension and compression and functionally graded properties of materials varying in a form of continuous function. For this purpose, we should resort to subareas model in tension-compression and layer-wise theory to simulate the problem studied, which inevitably complicates the analysis process further. The detailed steps for modeling and simulation are as follows.
(i) Establishment of entity structure The solid model of a FGPM cantilever beam is established, in which the length of the beam l is set to be 1 m, the section height h to be 0.2 m and the section width b to be 0.08 m.
(ii) Determination of neutral layer and tension-compression subarea Note that comparing to common piezoelectric cantilever beam, the subarea of tension and compression of beam under external loads is a newly introduced feature, we must determine the unknown neutral layer first, thus realizing the so-called tension-compression subarea. For this purpose, a set of functionally graded indexes should be chosen to determine the tensile height and the compressive one of the beams, according to the one-dimensional theoretical solution. For example, consider the case 1    As indicated above, the change of material properties as certain direction in ABAQUS cannot be defined as a continuous function, therefore according to conventional practice, we adopt layer-wise model to simulate the functionally graded properties of the materials. Without losing the computational accuracy, the beam is divided into a moderate number of layers; the physical parameters of the material in each layer are regarded as the same, thus indirectly realizing the continuous change of materials properties as the thickness direction if the layer numbers are enough. For this purpose, bounded by the neutral layer, the upper part and the lower part of the beam are divided equally, thus the beam is divided into 40 layers along the thickness direction, with each layer 5 mm thick, as shown in Figure 5. It is easy to see that there are 12 layers in the tensile zone and 28 layers in the compressive zone. Note that the coordinate origin now locates at the neutral layer.
continuous change of materials properties as the thickness direction if the layer numbers are enough. For this purpose, bounded by the neutral layer, the upper part and the lower part of the beam are divided equally, thus the beam is divided into 40 layers along the thickness direction, with each layer 5 mm thick, as shown in Figure 5. It is easy to see that there are 12 layers in the tensile zone and 28 layers in the compressive zone. Note that the coordinate origin now locates at the neutral layer.

(iii) Determination of properties of materials
The material constant at the neutral layer 0 z = is shown in Table 1, which may be input into ABAQUS directly. Material constants on other layers which do not locate at the neutral layer may be computed and input into the program, via the layering model of tension-compression subarea established in Step (ii).  (iv) Establishment of boundary conditions The left end of the beam is free and the right end is fully fixed, which agrees with the mechanical model shown in Figure 3. For this purpose, we need to define the fixed constrain on the right end of the cantilever beam in ABAQUS, including displacement and rotation.
(v) Mesh division In this simulation, an 8-node linear piezoelectric brick C3D8E is adopted, in which the size ratio is set to be 0.01 and the mesh size is 5 mm × 10 mm. C3D8E may well realize the simulation of a cantilever beam under the mechanical load and electrical load. (vi) Step module and adding loads An analysis step named Static General is established to apply the load. In this simulation, uniformly-distributed load form is considered only. For this purpose, a uniformly-distributed load (q = 1 N/m 2 ), along the negative direction of z axis, is applied on the upper surface of the beam.
Up to now, the modeling job has been finished. We note that in the above steps, the difference introduced by bimodular functionally graded properties embodies mainly in Steps (ii) and (iii), which makes the analysis more complicated. (vii) Operation and results output

(iii) Determination of properties of materials
The material constant at the neutral layer z = 0 is shown in Table 1, which may be input into ABAQUS directly. Material constants on other layers which do not locate at the neutral layer may be computed and input into the program, via the layering model of tension-compression subarea established in Step (ii).  The left end of the beam is free and the right end is fully fixed, which agrees with the mechanical model shown in Figure 3. For this purpose, we need to define the fixed constrain on the right end of the cantilever beam in ABAQUS, including displacement and rotation.
(v) Mesh division In this simulation, an 8-node linear piezoelectric brick C3D8E is adopted, in which the size ratio is set to be 0.01 and the mesh size is 5 mm × 10 mm. C3D8E may well realize the simulation of a cantilever beam under the mechanical load and electrical load.

(vi)
Step module and adding loads An analysis step named Static General is established to apply the load. In this simulation, uniformly-distributed load form is considered only. For this purpose, a uniformly-distributed load (q = 1 N/m 2 ), along the negative direction of z axis, is applied on the upper surface of the beam.
Up to now, the modeling job has been finished. We note that in the above steps, the difference introduced by bimodular functionally graded properties embodies mainly in Steps (ii) and (iii), which makes the analysis more complicated.

(vii) Operation and results output
After a job is established, submit the job and calculate and output the results. It should be noted here that although the numerical simulation is based on a three-dimensional case, the problem we study still attributes to two-dimensional plane problem. Therefore, we only output the results concerning plane problem, and meanwhile, it is also convenient to compare these outputs with our previous study [27] which is exactly the case of two-dimensional problem. Figure 6 shows the cloud diagram of the mechanical stresses, σ x , σ z and τ xz , the mechanical displacement, u and w, as well as the electrical displacement, D x and D z in a two-dimensional problem. After a job is established, submit the job and calculate and output the results. It should be noted here that although the numerical simulation is based on a three-dimensional case, the problem we study still attributes to two-dimensional plane problem. Therefore, we only output the results concerning plane problem, and meanwhile, it is also convenient to compare these outputs with our previous study [27] which is exactly the case of two-dimensional problem. Figure 6 shows the cloud diagram of the mechanical stresses, x σ , z σ and xz τ , the mechanical displacement, u and w , as well as the electrical displacement, x D and z D in a two-dimensional problem.

Comparison of One-Dimensional Solution and Two-Dimensional Simulation
Due to the fact that in the derivation of one-dimensional theoretical solution, some assumptions have to be introduced to obtain a simple but clear expression, the validity of one-dimensional solution should be further verified. For this purpose, we use the results from the two-dimensional numerical simulation to validate the rationality of the one-dimensional solution. Figure 7 shows the comparison results between one-dimensional theoretical solution and two-numerical numerical simulation. Figure 7a-c and Figure 7d-f show that the bending stress σ x and the electrical displacement D x at different cross sections x = 0.25l, 0.5l, 0.75l vary with the thickness direction, respectively; Figure 7g shows the deflection w(x) at z = 0. From Figure 7, it is easy to see that the two solutions curves are roughly close to each other, which indicates the validity of one-dimensional theoretical solution, to some extent. Besides, Tables 2 and 3 further show the relative errors of the one-dimensional solution and two-dimensional simulation concerning σ x at x = 0.25l and the vertical deflection w(x), respectively. It is obvious that the relative errors are within acceptable limits. Note that in Table 2, the value of ABAQUS simulation at the neutral layer z/h = 0 gives 0.2813 but zero (while the theoretical solution gives zero, according to Equation (11)). This error may be caused by the large discreteness of the finite element calculation itself, especially near the neutral layer. The neutral layer may be regarded as a computing-sensitive zone in which the tensile and compressive stresses change their positive or negative sign at this layer and for bimodular functional graded materials, coefficient of materials is continuous at this layer but their first-order derivative to z-direction is not continuous, see our previous study [29]. To obtain more accurate results, there is a need to increase the number of mesh division, especially at the neutral layer, thus enlarging the amount of computation.
It should be noted here, that since the one-dimensional theoretical solution is derived on a relatively simple case, the number of physical quantities obtained is relatively limited. Due to pure bending, for example, only σ x is derived while τ xz and σ z cannot be obtained in this way; similarly, only the deflection w(x) may be obtained in the one-dimensional solution while the so-called horizontal displacement u cannot be considered.  It should be noted here, that since the one-dimensional theoretical solution is derived on a relatively simple case, the number of physical quantities obtained is relatively limited. Due to pure bending, for example, only x σ is derived while xz τ and z σ cannot be obtained in this way; similarly, only the deflection ( ) w x may be obtained in the one-dimensional solution while the socalled horizontal displacement u cannot be considered.

Comparison of Two-Dimensional Numerical Simulation and Existing Solution
He et al. [27] derived a two-dimensional theoretical solution for FGPM cantilever beam with bimodular effect under combined loads. Therefore, it is interesting to compare the two-dimensional numerical results presented in this study with the existing solution from He et al. [27]. Figures 8-10 show the comparison results of the two-dimensional numerical simulation and the existing solution, in which the computational data of the existing solution are from [27]. Tables 4-6 further show the relative errors of the two-dimensional simulation presented in this study and existing work [27], in which σ x , τ xz , u, D x and D z at x = 0.25l are considered. The results show that the two solutions are approximately equal, and the relative error is also in the acceptable range, excluding a few points.

Comparison of Two-Dimensional Numerical Simulation and Existing Solution
He et al. [27] derived a two-dimensional theoretical solution for FGPM cantilever beam with bimodular effect under combined loads. Therefore, it is interesting to compare the two-dimensional numerical results presented in this study with the existing solution from He et al. [27]. Figures 8-10 show the comparison results of the two-dimensional numerical simulation and the existing solution, in which the computational data of the existing solution are from [27]. Tables 4-6 further show the relative errors of the two-dimensional simulation presented in this study and existing work [27], in which x show that the two solutions are approximately equal, and the relative error is also in the acceptable range, excluding a few points.   From the above results, it is easy to see that among the two-dimensional mechanical physical quantities, stresses σ x , τ xz and σ z as well as displacements u and w, their importance is different. For stress components, the bending stress σ x is still the dominant stress, the shear stress τ xz is secondary and σ z is almost negligible; for displacement components, it is obvious that only w is the interesting variables in our analysis while u is negligible small. This conclusion is consistent with the existing results. Similarly, among the two-dimensional electrical physical quantities, for example, D x and D z , the final result indicates D x > D z and in some cases even D x >> D z , this fact may well explain the rationality of the assumption D z ≈ 0 in the derivation of one-dimensional theoretical solution.  In addition, it should be noted here that, from Figures 9e and 10e, we may find that D z obtained by the two solutions have obvious difference, and for some points the maximum relative error is even greater than 50%. However, comparing with D x having the same attribute, D z is an unimportant quantity, as indicated above. Therefore, even if there is slightly big difference in D z , its influence on the whole problem is relatively limited, this is the reason why we obtained the one-dimensional theoretical solution by neglecting D z . Since the quantity itself can be ignored, the difference of this quantity is even more insignificant.

Comparison of One-Dimensional Theoretical Solution and Existing Solutions
In existing studies, Shi and Chen [17] gave a set of analytical solutions for the problem of a FGPM cantilever beam subjected to different loadings, and Yang and Liu [31] derived the solution for the problem of a piezoelectric cantilever beam under an end load. In this section, we will compare the one-dimensional solution derived in this paper with the two solutions mentioned above.
Before comparison, it is necessary to degrade these solutions into a one-dimensional case since they were obtained based on different conditions. For example, in [17], the effect of body force F z is considered, and only the elastic parameter S 33 changes along the z direction; while in this study, the different properties in tension and compression are taken into account, and all the material parameters change along the z direction. For the convenience of comparison, we let the body force F z from [17] be zero, and the elastic parameter S 33 from [17] be a constant, that is, where the meaning of these quantities above can be found in [17]. At the same time, we redefine some quantities from Section 3 as follows Thus, the problem solved in this study and the problem in [17] are now degraded to the problem of a piezoelectric cantilever beam under an end load, that is, the problem solved in [31]. The detailed degradation is shown as follows.
Firstly, let us degrade the solution obtained in [17]. Substituting Equation (46) into the Equations (14) and (17) in [17], we may obtain From Equation (48) and the Equation (11) in [17], it can be obtained that Thus, from Equations (48) and (49), the Equation (8) in [17] can be written as where " " represents the corresponding quantity derived from [17]. From Equations (48) and (49) as well as the Equation (13) in [17], we have When z = 0, the deflection equation at the axis of the piezoelectric cantilever beam can be obtained Next, let us degrade the one-dimensional solution derived in this study. Substituting Equation (47) into Equation (9), we have Since α 1 = α 2 = 0, we may obtain The integrations in Equation (20) may be simplified as From Equations (55) and (24), we obtain Integrating the above equation with respect to x will yield From Equation (26), we may obtain Thus Substituting Equation (56) into Equation (53), we finally obtain By comparing the above equation with Equation (50), it is easy to see that they are identical except the beam width b, which is considered to be unit 1 in [17]. In addition, from Equations (50) and (60), it is easy to find that σ x is the same as the corresponding expression, i.e., the Equation (29) in [31]. For the deflection w, Equation (59) in this study and the Equation (38) in [31] are identical, only with a slight difference in Equation (52). We note that s 11 (d 31 ) 2 /λ 33 for the piezoelectric materials PZT-4, as indicated in Section 3.1, thus s 11 − (d 31 ) 2 /λ 33 ≈ s 11 . From Equations (52) and (59), we may find, w ≈ w.

Evolution for One-Dimensional Theoretical Solution
For the one-dimensional solution of beams with different properties of materials, it is convenient to compare all kinds of solution due to the consistency in form. Among the solutions, we take two important physical quantities as our comparing objects, one is modulus of elasticity of materials, E * , and another is bending stiffness of beams, D * , since they are closely associated with final solutions. For example, if bending stiffness D * is known, the deflection is also easily determined by the classical Euler-Bernoulli equation in small-deflection case, 1 ρ = d 2 w dx 2 = M D * , where D * is associated with the modulus of elasticity E * of materials. For this purpose, Table 7 gives E * and D * , step by step, from classical beams to bimodular FGPM beams, via bimodular beam [33] and bimodular FGM beam [29]. Interestingly, this order may be called an evolution of material properties from classical beams to bimodular FGPM beams and, in turn, may also be called a regression from bimodular FGPM beams to classical beams. Table 7. Evolution for one-dimensional solution of beams (rectangular section).

Material Types of Beams Modulus of Elasticity Bending Stiffness
Classical beams E = const. 1 12 bh 3 E Bimodular beams [33] E + , Bimodular FGM beams [29] E + (z) = E 0 e α 1 z/h E − (z) = E 0 e α 2 z/h , where α 1 α 2 and E 0 is neutral layer modulus Bimodular FGPM beams (this study) It should be noted here that there is an important difference in A + 2 and A − 2 for bimodular FGM beams and bimodular FGPM beams. It is found that the integral function index has more than one negative sign for bimodular FGPM beams, i.e., e α i z/h in bimodular FGM beams is now changed as e −α i z/h in bimodular FGPM beams, this because, for bimodular FGPM beams, the original definition for material properties is based on flexibility coefficient s +/− 11 but stiffness coefficient E +/− , and they satisfy s +/− 11 = s 0 11 e α i z/h = 1/E +/− (z), thus generating the negative sign.

Discussion on Flexible FGPM Cantilever Beam
In many electromechanical devices, ultra large deflections of cantilevers are generally needed for the sake of application requirements. For example, Merupo et al. [34] investigated the flexoelectric response in soft polyurethane films and their use for large curvature sensing. More recently, Seveno and Guiffard [35] presented the realization of a cantilever-based PZT thin film deposited onto an ultra-thin aluminum foil as a substrate and showed that a very flexible actuator with low voltage-induced ultra large deflections can be obtained by this method. In the above studies, the large deformation analyses of piezoelectric structures are required, not only for static or dynamic problems but also for sensor or actuator models. Thus, via the one-dimensional solution obtained in this study, the deflection of flexible piezoelectric cantilever beam is easily obtained.
The bending stiffness of bimodular FGPM beams, D * , has been derived in the Section 3.1, which gives Here, if a large deflection bending is to be considered, the classical Euler-Bernoulli equation will be where ρ is still the curvature radius of the cantilever beam and w is the deflection. We note that here the curvature expression is mathematically precise and has not been the small curvature case mentioned above which reads 1/ρ = d 2 w/dx 2 . Via our previous studies concerning flexible cantilever beams made of classical materials [36], it is convenient to derive the deflection only by replacing the bending stiffness term D * in the solution obtained. This conclusion shows, from the side, the advantage of one-dimensional theoretical solution in predicting deformation of flexible piezoelectric cantilever structures.

Concluding Remarks
In this study, we used analytical and numerical methods to investigate a FGPM cantilever beams with different properties in tension and compression, in which one-dimensional theoretical solution was derived and two-dimensional numerical simulation was also performed. We made extensive comparisons to validate the rationality of the one-dimensional solution and two-dimensional numerical simulation in this study. The following main conclusions can be drawn.
(i) In the one-dimensional theoretical solution obtained, the presence of modulus of elasticity E * and bending stiffness D * for a bimodular FGPM beam may clearly describe the piezoelectrical effect on the classical problem without electromechanical coupling. The influences introduced by three important material coefficients, i.e., elastic coefficient, piezoelectric coefficient, and dielectric coefficient are all included in E * and D * .
(ii) In the two-dimensional numerical simulation, the layered model of tension-compression subarea opens the possibilities for the realization of numerical technique to the problem of FGPM cantilever beam with bimodular effect, although the software itself does not involve different properties in tension and compression of materials.