# Vibration Analysis of Piezoelectric Cantilever Beams with Bimodular Functionally-Graded Properties
## 4 Numerical Simulation
### Initial Modeling  
(i) Establishment of structures  
Suppose that the length \(l\) of a FGPM cantilever beam is set to be 50 mm, the section width \(b\) to be 10 mm and the section height \(h\) to be 2 mm.  
Determination of tension-compression subarea under initial displacement  
For the determination of the tensile and compressive height of the beam, a set of functionally-graded indexes should be chosen first in the light of our previous study on static analysis [28]. In this study, we consider the following two groups of gradient indexes: (a) \(\alpha_{1}=-2,\,\alpha_{2}=-3\) and (b) \(\alpha_{1}=2,\,\alpha_{2}=3\) to correspond to the two different cases: \(E^{+}(z)>E^{-}(z)\) as well as \(E^{-}(z)>E^{+}(z)\), respectively, as shown in Figure 3, in which \(E_{0}\) is the electric modulus value on the neutral layer, in case (a), \(h_{1}=0.6\) mm, \(h_{2}=1.4\) mm and in case (b), \(h_{1}=1.4\) mm, \(h_{2}=0.6\) mm.  
As is the case in most commercial software, there will be the limitations for implementing FGMs in ABAQUS; for example, it seems to be inconvenient to realize the change of material properties along a certain direction as a continuous function. To overcome the shortcomings, an alternative implementation was proposed in the context of the commercial finite element package ABAQUS [31]. In this study, however, the layer-wise model was still used to simulate the functionally-graded properties, since this practice is conventional and well-known. Without losing the computational accuracy, we divided the beam into a moderate number of layers along the thickness direction, the physical parameters of the material on each layer are considered to be the same, thus indirectly realizing the continuous variation of properties of materials along the thickness direction if the numbers for layering are sufficient. To this end, bound by the neutral layer, the upper and lower areas of the beam are equally divided into 40 layers along the thickness direction, each layer being 0.05 mm thick, as depicted in Figure 4. It is easy to see that in case (a), there are 28 layers in the compressive area and 12 layers in the tensile area, while in case (b), the layering is the opposite, that is, 12 layers in the compressive area and 28 layers in the tensile area. The coordinate origin is still on the neutral layer.  
(iii) Input of properties of materials  
The material PbZrTiO\({}_{3}\)-4 (generally abbreviated as PZT-4) is selected as our materials simulated. Table 1 shows the material constant on the neutral layer \(z=0\) which are directly input into ABAQUS.  
Figure 4: Sketch of layering on cross section of the beam (unit: mm): (a) \(\alpha_{1}=-2\), \(\alpha_{2}=-3\); (b) \(\alpha_{1}=2\), \(\alpha_{2}=3\).  
Material constants up or down the neutral layer may be computed and input into the program, by using the layering pattern of tension-compression established in Step (ii).  
Boundary conditions  
The left end of the beam is fully fixed and the right end is free; this agrees with the mechanical model presented in Figure 2. Besides, the upper and lower surfaces of the beam are open circuited.  
Mesh division  
An 8-node linear piezoelectric brick C3D8E is used, and the mesh size is set to be \(1\) mm\(\times\) 1 mm, in which the size ratio is 0.001 and also the global seed is set.