Static and natural frequency investigation of FGP beams considering thermal effects and design parameters

This study presents a static analysis and natural frequency analysis of functionally graded laminated piezoelectric beams based on the Euler–Bernoulli theory using the finite element method. A simple power law is used to vary all material properties across the thickness, except for Poisson’s ratio. The effect of laminate configuration and volume fraction index on the deflection and natural frequency of beams made of functionally graded piezoelectric materials (FGPM) is investigated, and the relationship between deflection and different volume fraction indices under thermal, electrical, and mechanical loads is explored. The study shows that there is a certain volume fraction index that maximizes or minimizes deflection. Additionally, the variation of natural frequency in relation to the power law index is examined. The findings of this research are useful for the development of sensors and actuators in different environments, and the appropriate operation point of the structure can be selected based on the behavior of the sensor or actuator of the beam.


Importance of the problem
In the past few decades, multilayer piezoelectric structures have become increasingly important in various fields of engineering as they can serve as sensors and actuators for monitoring the state or controlling vibrations of structures [1].However, such structures can also face several issues due to sudden changes in material properties at the junctions of individual laminates, such as residual stresses and delamination, especially in hightemperature environments [2].Stress concentrations between adjacent layers, creep at high temperatures, and failure at layer boundaries are also potential issues that can arise in multilayer panels subjected to mechanical, thermal, or electrical loading.A promising solution to address these limitations is the development of functionally graded piezoelectric materials (FGPMs), in which one or more layers have graded properties that smoothly vary through the structure, thereby reducing these shortcomings [3].

Background information of researches
The most commonly used piezoelectric sensors and actuators are laminated piezoelectric composites consisting of multiple layers, each layer having particular properties.Several studies have investigated the behavior of piezoelectric laminated composite structures.Robins and Reddy [4] investigated the static and dynamic interaction of the piezoelectric layer and the substructure layer of a clamped-free beam using Reddy's generalized laminated plate theory.Lee and Saravanos [5] developed a finite element formulation for a beam element with linear shape function.They discovered that there is a significant coupling between thermal and piezoelectric effect in the open circuit.Mukherjee and Chaudhuri [6] investigated the large deformation effect of piezoelectric composite materials.They compared the linear and nonlinear behavior of the carrier considering large deformations.Li et al. [7] studied the displacement of heterogeneous piezoelectric bimorph beam.They introduced a direct analytical solution to predict the displacement distribution in a heterogeneous piezoelectric bimorph beam.Varelis and Saravanos [8] analyzed the stability of piezoelectric shells under thermal loading.They derived formulation of constitutive equation of plates in curvilinear coordinate.The static bending, free vibration and dynamic response of FGPM actuators was analyzed by Yang and Xiang [9].The comparison between Euler-Bernoulli and Timoshenko beam theories illustrates that Timoshenko beam theory predicts large displacement.Xiang and Shi [10] investigated an FGP sandwich beam under thermal and electrical loading based on piezoelasticity theory.They found out that the electromechanical coupling must not be disregarded while the piezoelectric coefficient of the material is large.Behjat et al. [11] studied the static bending, free vibration and dynamic reaction of piezoelectric plates using the finite element method under mechanical and electrical loading.Non-linear behavior of FGPM plates were studied using FEM under various loadings by Behjat and Khoshravan [12].Lezgy et al. [13] introduced a new FE model for static, free vibration, and dynamic analysis of FGM piezoelectric beams using beam elements.Based on the strain gradient theory and Timoshenko beam theory, a size dependent beam model is developed by Li et al. [13].A model for simulating FGPM harvesters was introduced by Amini et al. [14].They found that the variation of power law index affects the sensory voltage and power of harvester.The three-dimensional Element-Free Galerkin method was developed by Mikaeeli and Behjat [15] in order to investigate the static behavior of thick FGPM rectangular plate with arbitrary material properties.They investigated the deformation, electric potential and electric field in the functionally graded piezoelectric plates.Nourmohamadi and Behjat [16] studied the static bending of FGPM plates under (thermo-electro-mechano) loading and different boundary conditions in order to find the specific volume fraction index that can be used in design of functionally graded piezoelectric plates.They discovered that there is a specific volume fraction index in which the plate deflection is maximized or minimized under thermal and electrical loading.Dai et al. [17] studied distribution of temperature, displacement, and stress of a rotating disk in the constant angular velocity under a coupled hygrothermal field by using a new numerical method.Bodaghi et al. [18] studied the effect of geometrical non-linearity and mechanical and thermal loadings on the active control of dynamic response of simply supported functionally graded beams under blast pulses using piezoelectric sensors and actuators.They found out that power law indices of FGP sensor and actuator can play a significant role.Bodaghi et al. [19] introduced the thermo-mechanical behavior of composite plates and used a new formulation.Kamarian and et al. [20] investigated free vibrations analysis of the FGSCNTR beams rested on Pasternak and studied the influences of four parameters in the CNT volume fraction relation on the natural frequencies.Salim et al. [21] investigated the thermal vibration of SMAHC cylindrical shells and found out that the orientation of SMAs in the layers plays an effective role in vibrations of SMAHC shells under thermal environments.Shen and Yang [22] investigated the large amplitude vibration of FG-FRC laminated beams and they found out that, temperature variation has a moderate effect on the natural frequencies of the hybrid laminated beam.Huang, Ding, and Chen [23] presented the analytical solutions for a functionally graded actuator subjected to constant electric potential for the two-dimensional plane stress problem.Pietrzakowski [24] discussed the comparison of changes in both the natural frequencies and resonant amplitudes depending on the PZT gradation.It was also seen that this actuator layers give a satisfactory operational effectiveness of the control system.Satyajit and Ray [25] investigated the geometrically nonlinear dynamic response of functionally graded laminated composite plates, where the numerical results showed that unlike the conventional laminated composite plates, the variations of stresses across the thickness of the FG laminated composite plates are smooth and continuous.Wu and Liu [26] introduced a semi-analytical method for three dimensional analyses of composite FGPM plates and shells.Ridha et al. [27], employed a higher order shear deformation beam formulation having three variables without using of shear correction factor in order to analyze forced vibrational of shear deformable functionally graded (FG) nanobeam under partial dynamical load, consequently they found out as the dynamical force moves away from the beam edges, the dynamical deflections increase.
Ridha et al. [28] investigated nonlinear thermal stability behaviors of elastic nanobeams with piezo-magnetic properties.They employed Nonlocal theory for mathematical formulating and they found that the buckling temperature reduced via higher rates by increase of piezoelectric constituent volume.Jahanghiry et al. [29] studied the mechanical behavior of an FGM microgripper under the effect of DC voltage and temperature variation the found that the static pull-in voltage of the FGM microgripper, which is initially subjected to the initial temperature variation, decreases by increasing the ceramic volume, and instability occurs at lower DC voltage.Maleki and Mohammadi [30] investigated study effect of piezoelectric patches and crack parameters on stability of the cracked FGM columns their Results show that the crack significantly reduces buckling load of the column.

Description of current problem
This article investigates the thermo-electro-mechanical analysis of the FGPM beam using FEM while the FGPM beams features are considered and investigated in most of the researches, in present paper, a new type of multilayer piezoelectric beam introduced by adding two piezoelectric layers to FGPM beam.The introduced beam was investigated under variety of loading and boundary conditions.There are noticeable papers dealing with the FGPM beams, however majority of them overlooked the impact of different volume fraction index on the static deflection in thermal and electrical loadings.In this paper the optimum response of multi-layered FGP beams studied and the results shows that there is a particular power law index in which the deflection of beam maximized or minimized so that be used for design of sensors and actuators.Although the natural frequency of the beam is affected by altering the power law index so the effects of volume fraction index on the deflection of the FGPM beam become important when the FGPM device has to perform efficiently in electrical or thermal environments, the results of this study could certainly take steps to address these issues.

Methods
In this paper, the static deflection and natural frequency of multilayer FGPM beams were investigated based on the Euler-Bernoulli beam theory.The deflection of the beam was calculated under thermo-electro-mechanical loading.The FE method based on linear elements was used to study the beam.The beam equation was obtained using the equilibrium equation of mechanical and electrical imbalance forces.The beam behavior was studied for different laminate configurations and volume fraction indices, and an optimal volume fraction index was determined to design multilayer smart structures in a thermal environment.The first three natural frequencies of the multilayer FGPM beam were determined for different boundary conditions and different laminate configurations as a function of the volume fraction index.Our aim is to find a specific value of volume fraction index that the beam deflection is maximized or minimized.

Functionally graded piezoelectric beams
Various models have been presented to simulate material distribution throughout the thickness of the beam.However, usually the power law distribution is used to simulate the distribution of material properties.For FGM structures made of two different materials, material "a" on the top side and material "b" on the bottom side of the beam, the effective properties of the FGPM profile in the thickness direction of the beam can be described as follows [31]: The effective material property: of the FGPM beam, denoted as "P eff " and "P a " is the upper surface property of the FGPM beam, "P b " is the lower surface property of the FGPM beam, and "V a " is the volume fraction of the functionally graded piezoelectric upper material, which can be expressed as The volume fraction of the bottom surface of the beam, denoted as V b and "n" is the power law index.Thus, the material properties can be written as (1) The distance from the neutral axis of the beam, denoted as "z" and "h" represents the beam thickness.In this article, all of the material properties differ throughout the thickness by following the simple power law distribution (Eq.4) excluding Poisson's ratio.

Euler-Bernoulli beam theory
The displacement field for Euler-Bernoulli beam is [32] Here, u and w are the displacements in the x and z directions, respectively.The dis- placement field in vector form can be written as follows: In Eq. ( 6), the variable ℜ defined as follows: Additionally, the strain field can be written as

Force and moment resultants formulation
In each lamina, the stresses can be integrated throughout the laminate thickness to determine the resulting forces and moments.Let 'n' represent the number of plies in the laminate.Therefore, suppose a laminate consists of 'n' plies and each ply has a thickness of "tk".The resultant force and moment in the laminate are written in matrix form as [33]: The normal force per unit length, defined as N x and M x is the bending moment per unit length.
The global stress in each lamina is (4)

Constitutive equations of piezoelectric material
Based on the linear theory of piezoelectricity, the constitutive equation for piezoelectric beam can be written as [34] The stress and strain tensors, denoted as σ ij and ε kl , respectively and C E ijkl is the elasticity matrix, e ijk is the piezoelectric coefficient tensor, E k is the electrical field vector, ij is the thermal expansion tensor, �θ is the temperature difference from the reference tem- perature, D l is the electrical displacement vector, ξ lk is the electrical permittivity, and p l is the pyroelectric vector.Also, ij can be written as In Eq. ( 12), α kl is the thermal expansion coefficient.The thermal load considered in this analysis is assumed to be constant at the top and bottom surfaces.The temperature field is changed only in the thickness direction.Therefore, the heat transfer equation is [35] The thermal conductivity, denoted as "K" that is assumed to obey the simple power law distribution.The temperature field is obtained throughout the thickness of the beam by solving Eq. ( 13) as follows:

Weak form of equilibrium equations
Using the divergence theorem and neglecting the inertial effects, the variational figure of the equations of motion is written as follows [36]: The vectors ψ u and ψ e shows unbalances between internal and external forces.Also, parameters b, τ , q and v describe the body force, surface traction on the bounding sur- face and φ is the electric potential.(10) The combination of Eq. 15 and Eq.11 the final equation can be achieved as Here N describes the number of piezoelectric layers, G is electrical permittivity.In addi- tion, T is the pyroelectric component.The other matrices can be represented as

Finite element formulation
In this article, two-node linear element is applied in order to analyze the FGPM beam.Each node has four degrees of freedom { � u , φ } .The displacement field represented as: The Lagrangian shape function " ψ " is used to interpolate "u 0 ", along with Hermit cubic shape function "N" which is interpolating function for deflection "w".These approximations are defined as ( 16) And for each element, the displacement vector is interpreted as In this case, { u e } and [N u ] are the nodal displacement and interpolation matrixes that are written as The strain relation could be written as Here L and B u are defined as [R a ] and [R b ] represent the strain operators.These operators can be written as Using the linear Lagrangian interpolation function N φ , the potential φ and the elec- tric vector E are interpreted the nodal variables as [37] And the electrical field can be written as The electric field varies in the direction of thickness as (E z = ∂ϕ ∂z ) .Assuming that the electric potential remains constant across the beam thickness, the electric field can be defined as follows: (19) In this case, thickness of piezoelectric layer is shown by "h p ". Placing Eq. 16 in Eq. 23 and integrating throughout the thickness and length direction of the beam, the equation of motion of FGPM beam in matrix form can be acquired [5] Here submatrices K ee , K de , and K dd demonstrate the permittivity stiffness, piezoelectric and elasticity matrices, respectively.
The sensory electrical potential is obtained as By combining Eq. 31 and Eq. 30, the following equation for displacement could be obtained: Here [F ] is the mechanical force vector.In addition, thermal, coupled thermal, and electri- cal load vectors can be defined as The abstract form the of final equation of motion of FGPM beam can be written as Equation 31 can be written as By solving above equations, displacement field in the FGPM beam can be obtained.In addition, by solving the Eigen-value problem, the natural frequencies of the beam are calculated: (29)

Comparison studies
To substantiate the results obtained in the present study, four examples are used to compare the data with the results of Lee and Saravanos [5], Yang and Xiang [9], and Phung et al. [38]: Example 1: Multiple piezoelectric layers consider the present example under a linear thermal gradient (120 °C at the bottom and 20 °C at the top) acting on the beam.Then, two different electrical voltages (100 and 200 V) are applied to the beam.Although the beam has a symmetrical laminate configuration including the piezoelectric layers, the applied temperature gradient causes the beam to bend.The beam has a length of (L = 25.4 cm) and a width of (b = 2.54 cm) and consists of different orientations of graphite/epoxy layers connected to piezoelectric lay-  ers, with each graphite/epoxy layer and piezoelectric layer having a thickness of (hl = 0.0127 cm).The deflection of the beam is given as 10 (w/b) without dimensions.From Fig. 1, it can be seen that the present analysis agrees very well with the results of Lee and Saravanos [5].
Example 2: In this section, a clamped-free FGPM beam is considered.The top of the beam is made of PZT-4 and the bottom is made of PZT-5.The material properties of the piezoelectric materials are based on a simple power law.Table 1 shows the deflection of the tip of the FGPM beam under mechanical load q = 10KN/m2 and electrical load V = 20 V in the volume fraction index of n = 0.2.The non-dimensionalized results are shown as (w/h).It is observed that the results obtained are in good accordance with the data reported by Yang and Xiang [9].Example 3: To verify the results, the data are compared with the results of Yang and Xiang [9].An FGPM beam consisting of (PZT-4 and PZT-5) is assumed.The natural frequencies of the beam are determined and compared with the results reported by Yang and Xiang.Table 2 shows the results for the natural frequency in this work and in Yang and Xiang.It can be seen that these results agree well with those of Yang and Xiang.
Example 4: A piezoelectric bimorph clamped-free beam with a width of 5 mm, a length of 100 mm, and a thickness of 1 mm is considered.The material of the piezoelectric beam consists of two PVDF layers.The beam is loaded with different voltages (50 V and 200 V).The deflection of the beam is shown in Fig. 2.These results show high match with the data of Phung et al. [38].Table 3 shows the errors reduction by increasing the number of element, thus in order to convergence the results and demonstrate the independence from the mesh, the results compared to the deflection of the tip of the FGPM beam under mechanical load q = 10KN/m 2 and electrical load V = 20 V in the volume fraction index of 0.2 which mentioned before in Table 1.Consequently, in all the calculation in the present study 25 elements are used.

Static analysis of multilayered FGPM beams
In this section, the static bending of functionally graded piezoelectric beam under electrical, thermal and mechanical loading is studied.The length to thickness ratio of the beam is L/h = 25.In this beam the thickness of each piezoelectric layer is 5% of total thickness of the beam "h" and the thickness of FGPM layer is 90% of total thickness of the beam.Figure 3 shows the laminate configuration of such beam.The material properties of mid layer obey the power law distribution, in which the bottom surface is made up of PZT-5 and the top surface is made up of PZT-4.The material properties of piezoelectric materials are listed in Table 4.In which Q 11 and Q 55 is the stiffness of the beam and d ij is the electrical displacement and "p' is the pyroelectric coefficients.The beam deflection is defined as, w = w h , where "h" and "w" are the beam thickness and the beam deflection, respectively.In this paper, for simplification, the three types of laminate configuration of the beam are considered and shown by "A" for FGPM beam with PZT-4 and PZT-5 layers, "B" for FGPM beam with both PZT-5 layers and "C" for FGPM beam with both PZT-4 layers.

Multilayered FGPM beam under mechanical loading
Figure 4a, b illustrates the dimensionless central deflection of a simply supported and tip deflection of a clamped-free Multilayered FGPM beam under uniform mechanical loading (q = 5KN/m) for all type of laminate configuration.According to the observation by raising the value of volume fraction index, the beam deflection is increasing.This phenomenon can be explained as the elastic modulus of PZT-5 piezoelectric is less than the elastic modulus of PZT-4.In addition, the deflection in the beam, type "B" is more than the other types of configurations.

Multilayered FGPM beams under uniform thermal loading
In this section, the beam in both bottom and top sides undergo a uniform (ΔT = 10 °C) thermal loading.Figure 5a, b describes the maximum dimensionless deflection of the simply supported and dimensionless tip deflection of multilayered FGPM beam versus volume fraction index in different types of configurations, respectively.From Fig. 5a, b, it is deduced that by increasing the volume fraction index, in types "B" and "C", the deflection is increasing initially and then decreasing.It is seen that there is a particular value of volume faction index in which the beam deflection is maximized in both boundary conditions "n = 1.5".In type "B" for both boundary conditions, the deflection is decreasing at the first, until the deflection is getting to zero in "n = 0.126", and from this point the deflection is increasing and getting maximized in "n = 1.5", and then, is decreasing.This phenomenon can be justified by the diversity in the value of the thermal expansion coefficient of PZT-4 and PZT-5 materials that can cause this type of deflection in different volume fraction indices.

Multilayered FGPM beams under thermal gradient loading
In this section the beam undergoes a 10 °C thermal gradient load, in which the bottom surface temperature is 10 °C and the top surface temperature is 0 °C.Figure 6a, b depicts the maximum deflection of FGPM beam that has simply supported or clamped boundary condition.From Fig. 6a, b, it can be induced that for both boundary conditions by increasing the power law index, the deflection is reducing at the first and then, increasing.There is a particular value of volume fraction index that the deflection is minimized.The value of this power law index in simply supported boundary condition is "n = 0.25" for type "B" and "C", and n = 0.15" for type "A".In case of clamped-free boundary condition, the value of this index is "n = 0.75″ for type "B" and "C", and n = 0.6″ for type "A".

Multilayered FGPM beams under electrical loading
In this section, an electrical loading (V = 100 V) is applied on the beam.Figure 7a, b shows the maximum dimensionless deflections of the multilayered FGPM beam with simply supported and clamped-free boundary condition in comparison with volume fraction index for different types of configuration, respectively.From Fig. 7a, b, it is inferred that by growing the power law index, deflection is raising first, and then, reducing.Also, there is a particular value of volume faction index "n = 1.6" that the deflection is maximized.It can be explained by the difference in the value of the piezoelectric constants of the constituents (PZT-4 and PZT-5) that can cause this type of deflection.

Beams under electrical and thermal loading
In the this section, The uniform thermal loading (ΔT = 10 °C) and electical loading (V = 100 V) is applied to the multilayered FGPM beam.In Fig. 8a, b, it can be seen that for both boundary conditions, the deflection is decreasing at the first until it reaches zero around "n = 0.30" for type "A" and "n = 0.1" for type "B" and "C" of laminate configuration, respectively and then, the deflection is increasing until it reaches its maximum value and then it is decreasing.The deflection is maximized around "n = 0" for model "A" and "n = 1.6" for model "B" and "C" of laminate configuration.
As it is seen from Figs. 7 and 8, an optimum point of deflection is occurred due to thermal and electrical loading.This phenomenon is deeply depend on the Young modulus, thermal expansion and piezoelectric coefficients of the materials of the beam.It means that by increasing of power law index, the thermal expansion coefficients decreases while Young modulus increases and this event lead to optimum point

Natural frequency analysis of multilayered FGPM beams
In this section, the natural frequency of multilayered FGPM beam is obtained for different boundary conditions, laminate configuration and volume fraction indexes.The results are non-dimensionalized as The first three natural frequencies of multilayered FGPM beam for different laminate configuration and bounady conditions is obtained and presented in Tables 5, 6, and 7 for  type "A", type "B", and type "C" respectively.It is observed that by increasing the value of volume fraction index, the natural frequency is decreasing in the multilayered FGPM beam and this phenomena can be explained by the variation of modulus of elastisity by increasing of power law index.Figures 9, 10, and 11 demonstrate the first natural frequency of the FGPM beam in different bounday conditions and various laminate configuration.It is observed that for differenet boundary conditions the type "C" of laminate configuration that consists of FGPM beam with two PZT-4 layers has higher natural frequency and the type "B"

Dimensionless frequency
Clamped-Free Boundary condition Type B Type C Type A Fig. 9 The first dimensionless natural frequency for clamped-free boundary condition laminate configuration has lower natural frequency.The reason is related to the modulus of elasticity of constituent.

Conclusions
In this paper, a new configuration of multilayered FGPM beam was introduced and also the static and natural frequency of multilayered FGPM beam was investigated in thermoelectro-mechanical loading and different laminate configurations.The material properties are assumed to obey the simple power law distribution.The static deflection of the beam was calculated by using finite element method based on two-node linear element and the natural frequency was obtained by solving the related Eigen-value problem.The multilayered FGPM beams behavior was studied under different loading type and volume fraction indices.It was inferred that there is an spacific value of volume fraction Fig. 11 The first dimensionless natural frequency for simply supported boundary condition index that the beam deflection is maximized or minimized.For a multilayered FGPM beam under electrical or uniform thermal loading, there is an specific volume fraction in which the deflection is maximized.Also, for an FGPM beam under thermal gradiant loading, there is an spacific volume fraction in which the deflection is minimized.This feature can be used as a design criterion in the selection of appropriate power law index for the FGPM sensors and actuators based on the maximum or minimum deflection.

Fig. 1
Fig. 1 Non-dimensional deflection of the beam under thermal gradient loading and different applied voltage 0, 100, and 200 V

Fig. 2
Fig. 2 Deflection of the piezoelectric bimorph beam under different applied voltages

Fig. 3
Fig. 3 Illustration of problem schematic and three types of laminate configuration

Fig. 4 a
Fig. 4 a Dimensionalized maximum deflection of center of the multilayered FGPM beam under mechanical loading q = 5KN/m in simply supported boundary condition.b Dimensionalized maximum deflection of free end of a clamped-free multilayered FGPM beam under mechanical loading q = 5KN/m

Fig. 5 a
Fig. 5 a Dimensionalized maximum deflection of center of the multilayered FGPM beam under uniform thermal loading for different types of configuration versus power law index in simply supported boundary condition.b Dimensionalized maximum deflection of free end of clamped-free multilayered FGPM beam under uniform thermal load for different types of configurations versus power law index

Fig. 6 a
Fig. 6 a Dimensionalized maximum deflection of center of the multilayered FGPM beam under thermal gradient loading for different types of configuration versus power law index in simply supported boundary condition.b Dimensionalized maximum deflection of free end of clamped-free multilayered FGPM beam under thermal gradient loading for different types of configurations versus power law index

Fig. 7 a
Fig. 7 a Dimensionalized maximum deflection of center of simply supported multilayered FGPM beam under electrical loading V = 100 V. b dimensionalized maximum deflection of free end of clamped-free multilayered FGPM beam under electrical loading V = 100 V

Fig. 8 a
Fig. 8 a Maximum dimensionless central deflection of simply supported multilayered FGPM beam under uniform thermal loading (ΔT = 10 C) and electrical loading (V = 100 V). b Maximum dimensionless tip deflection of clamped-free multilayered FGPM beam under uniform thermal loading (ΔT = 10 C) and electrical loading (V = 100 V)

Fig. 10
Fig.10The first dimensionless natural frequency for clamped-clamped boundary condition

Table 1
Maximum dimensionless deflection of free end of a clamped-free FGPM beam under electrical and mechanical loading

Table 3
Calculation errors in different number of elements

Table 5
The first three dimensionless natural frequencies for the FGPM beam type "A"

Table 6
The first three dimensionless natural frequencies for the FGPM beam type "B"

Table 7
The first three dimensionless natural frequencies for the FGPM beam type "C"