# A Unified Analysis for the Free Vibration of the Sandwich Piezoelectric Laminated Beam with General Boundary Conditions under the Thermal Environment

Authors: Guohua Gao; Ningze Sun; Dong Shao; Yongqiang Tao; Wei Wu
Journal: Shock and Vibration
 (2021)

---

## Abstract

This article mainly analyzes the free vibration characteristic of the sandwich piezoelectric beam under elastic boundary conditions and thermal environment. According to the first-order shear deformation theory and Hamilton’s principle, the thermo-electro-elastic coupling equations of the sandwich piezoelectric beam are obtained. Meanwhile, elastic boundary conditions composed of an array of springs are introduced, and the displacement variables and external potential energy of the beam are expressed as wave functions. By using the method of reverberation-ray matrix to integrate and solve the governing equations, a search algorithm based on golden-section search is introduced to calculate the required frequency parameters. A series of numerical results are compared with those reported in literature studies and obtained by simulation software to verify the correctness and versatility of the search algorithm. In addition, three parametric research cases are proposed to investigate the frequency parameters of sandwich piezoelectric beams with elastic restraint conditions, material parameters, thickness ratio, different temperature rises, and external electric potential.


---

## Body

## 1. Introduction

Piezoelectric ceramics are an electronic ceramic material with positive and negative piezoelectric effects; they have been widely used in sensors, energy harvesting devices, and smart actuators. Given their piezoelectricity, dielectric nature, elasticity, and other characteristics, piezoelectric ceramics are highly transformable in terms of size, external parameters, and lamination methods. Piezoelectric materials and general materials are usually laminated or pasted in engineering to achieve improvement in their performance. The ultrasonic motor widely used in the aerospace and automobile field is composed of two layers of piezoelectric sheets and a base metal material in the middle layer. The structure under such high-temperature and pressure environment is prone to overloading and damage. Therefore, a comprehensive understanding for vibration characteristics of the piezoelectric laminated structure in a thermal environment is necessary to reduce resonance through the rational design of structural parameters.Given that the aspect ratio of the piezoelectric laminated structure is generally greater than 3.5, it can be suitably simplified into a beam model. The related beam theories mainly include elasticity theory [1, 2], Euler–Bernoulli theory which is also named classical beam theory (CBT) [3, 4], first-order shear deformation theory (FSDT) [5–8], and high-order shear theory [9, 10]. The early CBT was suitable for studying the characteristics of slender beams. Given that the theory disregards the influence of transverse shear deformation, it is unsuitable for studying medium-thickness beams. As a remedial measure, the Timoshenko beam theory was obtained by adding shear deformation to the Euler–Bernoulli equation. Thus, the rotation of the cross section is caused by deflection and transverse shear. Different cross sections are needed to introduce the shear correction factor to reduce errors. For the higher-order shear theory without using the shear correction factor, different shape functions are introduced for the shear strain, resulting in the more accurate calculation of laminated plate problems.Given the various beam theories mentioned above, numerous researchers have achieved remarkable results in the analysis of piezoelectric beams. From the perspective of their shape, beams can be divided into bimorph beams [11–13], functionally graded piezoelectric beams (FGPM) and sandwich beams [14–22], nano-microbeams [23–27], and curved beams [28, 29]. The above literature shows that the transition from the traditional composite material lamination theory to the piezoelectric material lamination structure has been relatively complete. Carrera et al. [30] explained in detail the derivation process of piezoelectric and governing equations when smart materials are applied to beam structures and perfected the constitutive equation of the piezoelectric laminated structure under the above deformation theories. On this basis, the analysis of the dynamic characteristics of piezoelectric beams in various coupling environments has become an important research topic in recent years. Annigeri et al. [31] used the in-plate finite element method (FEM) to analyze the vibration characteristics of magneto-electroelastic beams based on CBT. He et al. [32] studied the free vibration of thermo-magneto-electro-nanobeams based on the wave solution method and carried out a series of parametric research including slenderness ratio, thickness, temperature rise, and external potential under different boundary conditions. Yang and Xiang [33] used the differential quadrature method (DQM) to study the FGPM beam under thermoelectric coupling. Lee [34] studied the displacement and stress responses of the FGPM beam with layer-wise finite element formulation and analyzed the effects of different loads on the deflection and stress of bimorph actuators. Benjeddou et al. [35] used the strain gradient theory to study the dynamic responses of the cantilever beam and discussed the effect of structure size and electric potential. Summarizing the results of these researchers, the external conditions for piezoelectric beams generally include the thermal environment, electric field, and magnetic field. However, these studies mostly focused on traditional functionally graded structures or microstructures, and a relatively limited number of research studies centered on sandwich structures under the above conditions. Therefore, the dynamic characteristics of sandwich piezoelectric beams under a coupling environment must be analyzed, and a parameterized study on their geometric properties and external conditions must be conducted.The analysis of the free vibration of this type of structure can ultimately be attributed to the problem of solving differential equations. In addition to various beam theories mentioned above, various analysis methods are used to solve the dynamic characteristics of composite material laminated beams under coupling fields; such methods include the FEM [36–40] and DQM [41–43], modal superposition method [44, 45], transfer matrix method [46–48], Galerkin method [49], Ritz method [50], and pseudo-arclength continuation method [51]. However, the modal superposition method has the problem of a large amount of calculation and low accuracy, the transfer matrix method also has the problem of insufficient high-order frequency accuracy, and the differential quadrature method is not convenient when calculating general boundary conditions. The research on structural boundary conditions also stays in traditional clamped, simply supported, and free. In actual engineering, there are more elastic connection ways such as bolt, rivet, and glue connections instead of the above methods. Therefore, an artificial boundary realized by changing the spring stiffness coefficient is introduced to obtain more accurate calculations. In addition, an analysis technique based on the frequency domain named the method of reverberation-ray matrix (MRRM) was proposed by Howard and Pao [52], which was used to solve the problem of elastic wave propagation in frame structures. The basic idea is to establish a pair of coordinate systems on the structure, derive the phase and scattering relationships of the structure through the dual relationship under two coordinate systems, and obtain the overall MRRM equation containing a series of dynamic properties of the structure. Guo et al. [53] discussed the MRRM in the functionally graded layered medium formula and equation-solving method. Miao et al. [54] used the MRRM to establish the frequency response function of the displacement in the frequency domain. The frequency is determined by the frequency response curve, and the mode shape can be calculated by the joint matrix from the system control equation. Thus, the MRRM is highly suitable for solving the dynamic characteristics of related structures. The scattering matrix is an all-constant matrix resulting from the dual coordinate system, its numerical value is relatively stable, and the high-order frequency can be accurately solved. However, limited research studied the dynamic characteristics of sandwich piezoelectric structures with arbitrary boundary conditions based on the MRRM.As in previous summaries, this article uses the MRRM to analyze the free vibration characteristics of the sandwich piezoelectric beam with elastic boundary conditions and coupling field. In accordance with the FSDT and Hamilton’s principle, the thermo-electric-elastic coupling equation and the boundary coordination relationship are obtained. Three sets of mechanical springs are used to implement any combination of the boundary conditions. On these bases, the frequency results are compared with the reported literature studies and software. For enriched practicability and flexibility of this theory, a series of parameters that affect the vibration characteristics, including the elastic modulus, thickness ratio, ply orientation of the middle metal layer, elastic boundary conditions, temperature rise, and external electric potential, is investigated.


## 2. Mathematical Model

### 2.1. Model Description

Figure1 shows a sandwich piezoelectric beam with piezoelectric sheets as the top and bottom layers and a general metal material layer in the middle. The length and thickness are set as L and h, respectively. The piezoelectric and metal layers have a thickness of hp and hm, respectively. Each layer is assumed to be perfectly pasted, and the influence of the adhesive layer is disregarded. For simulation of the boundary force at both ends of the beam, there are three sets of mechanical springs, including two sets of displacement springs and a set of rotary springs, which are set in each group. Furthermore, a set of electric displacement constraints is considered. In the following derivation, the piezoelectric material equation will be described first. For general metal layers, only the piezoelectric term in the equation needs to be removed. Moreover, zk + 1 and zk indicate the distance from the top and bottom of the kth layer to the middle face of the beam, respectively.Figure 1
Piezoelectric composite laminated beam with elastic boundary.
(a)(b)


### 2.2. Kinematic Equation

The external electric potentialϕ and temperature rise ΔT are introduced to describe the thermo-electro-elastic coupling equation of the sandwich piezoelectric beam. According to the assumption of the FSDT, the displacement of any point can be written as follows:(1)uxx,z,t=u0x,t+zθx,t,uzx,z,t=w0x,t.In (1), u0,w0,θ are the displacements and rotation angle of a node in the middle plane. The variables of external electric potential are given by the following [32]:(2)Φ¯x,z,t=−cosαzΦx,t+2zΦ0h,where α=π/h and Φ0 is the initial external electric potential.According to (1) and (2), the mechanical strains and electric field are, respectively, obtained as follows:(3)εxx=∂u∂x,=∂u0∂x+z∂θ∂x,γxz=∂w∂x+∂u∂z,=∂w0∂x+θ,Ex=−∂Φ¯∂x,=cosαz∂Φ∂x,Ez=−∂Φ¯∂z,=−αsinαzΦ−2Φ0h.Based on the stress-strain relationship of orthotropic materials and thermal stress, the constitutive equations of thekth layer from the sandwich piezoelectric beam can be written as(4)σxxk=c¯11kεxx−e¯31kEz−c¯11kα1kΔT,τxzk=c¯55kγxz−e¯15kEx,Dxk=e¯15kγxz+s¯11kEx,Dzk=e¯31kεxx+s¯33kEz+p¯3kΔT,where the overline indicates the transformed coefficient. c¯ij is the elastic constants; e¯ij and s¯ij are piezoelectric and dielectric constants, respectively. Notably, e¯ij and s¯ij are equal to zero when the kth layer is a general metal material; these constants can be specifically expressed as follows [23]:(5)c¯11=c11−c132c33,c¯55=c55,s¯11=s11,s¯33=s33+e332c33,e¯15=e15,e¯31=e31−e33c13c33,p¯3=p3−e33p3c33,where σxx and τxy represent the normal and shear stresses, respectively. Dx and Dz represent the electric displacement. Dy and Ey are neglected, and σyy=σzz=τxy=τyz=0. α1 and p3 correspond to the thermal expansion and pyroelectric constants, respectively. Thus, the thermal strain can be written as(6)εxxTk=α¯1kΔT.By integrating the mechanical stress in (4) in the thickness direction, the resultant expression of stretching stress, moment, and transverse force can be calculated with the following:(7)Nxx=∑k=1n∫zkzk+1σxxkdz,Mxx=∑k=1n∫zkzk+1zσxxkdz,Qxz=∑k=1n∫zkzk+1τxzkdz.Similarly, the corresponding stiffness coefficients can be obtained in the same integration way:(8)A11,B11,D11=∑k=1n∫zkzk+1c¯11k1,z,z2dz,A55=∑k=1n∫zkzk+1c¯55kdz,E31,F31=∑k=1n∫zkzk+1e¯31k1,zαsinαzdz,E15=∑k=1n∫zkzk+1e¯15kcosαzdz,R11=∑k=1n∫zkzk+1s¯11kcos2αzdzc,R33=∑k=1n∫zkzk+1s¯33kα2sin2αzdz,where κ represents the shear correction factor that is determined by a cross-sectional shape. Aij, Bij, and Dij refer to the mechanical stiffnesses, and Eij and Fij denote the electromechanical stiffnesses. Rij stands for the dielectric stiffness [30]. Therefore, we can obtain the constitutive equations of the piezoelectric laminated beam as follows:(9)Nxx=A11∂u0∂x+B11∂θ∂x+E31Φ,Mxx=B11∂u0∂x+D11∂θ∂x+F31Φ,Qxz=κA55∂w0∂x+θ−κE15∂Φ∂x,∫−h/2h/2Dxcosαz=E15∂w0∂x+θ+R11∂Φ∂x,∫−h/2h/2Dzαsinαz=E31∂u0∂x+F31∂θ∂x−R33Φ.


### 2.3. Governing Equations and Geometric Boundary Conditions

The coupling equation under thermal and electric fields of the sandwich piezoelectric beam is obtained by Hamilton’s principle defined by the following:(10)∫t1t2∫VδUM+UT+UE−KdVdt=0,where δ represents the variation and UM indicates the virtual strain energy and virtual electric internal energy calculated by the following [30]:(11)UM=∑k=1n∫0L∫zkzk+1σxxMεxxM+τxzγxz−DxEx−DzEzdzdx.Substituting equations (3), (4), and (7) into (11) and performing variational operations on the above equation, the expression of total virtual internal work is shown as(12)δUM=∫0LNxx∂δu0∂x+Mxx∂δθ∂x+Qxz∂δw0∂x+δθ∫−h/2h/2−Dxcosαz∂δΦ∂x−DzαsinαzδΦdzdx.UT and UE represent the external thermal potential energy and electric potential energy, respectively. In this article, only the axial strain component should be considered, and it can be expressed as(13)εxxT=12∂w∂x2=12∂w0∂x2.Thus, the external thermal and electric potential energies can be written with the following equations:(14)δUT=−∑k=1nm∫0L∫zkzk+1σxxTkδεxxTkdzdx=−∫0LNxxT∂w0∂xδ∂w0∂xdx,δUE=−∑k=1np∫0L∫zkzk+1DzkδεxxTkdzdx=−∫0LNxxE∂w0∂xδ∂w0∂xdx,where NxxT and NxxE are the normal forces as(15)NxxT=c11kεxxTk,NxxE=−2e¯31ϕ0.From the beam displacement field in (1), K represents the kinetic energy of the sandwich piezoelectric beam:(16)∫t1t2δK=−δ∫t1t2∫0L∫−h/2h/2ρku˙2+w˙2dzdx,where the material density of the kth layer is expressed as pk, and the superscript indicates the integral over time.(17)δK=−∑k=1n∫0L∫zkzk+1ρku¨δu+w¨δwdzdx=∑k=1n∫0Lρku˙02+zθ˙2+2zu˙0θ˙+w˙02dx=∫0LI0u¨0δu0+I0w¨0δw0+I1θ¨δu0+I1u¨0δθ+I2θ¨δθdx,where I0,I1,I2 are the kth layer’s inertia coefficients defined by(18)I0,I1,I2=∑k=1n∫zkzk+1ρk1,z,z2dz.Equations (12), (14), and (17) are substituted into equation (10) to derive the governing equations of the sandwich piezoelectric beam under the coupling environment as follows:(19)0=−∫0L∂Nxx∂x−I0u¨0−I1θ¨δu0∂Qxz∂x−NxxT+NxxE∂2w0∂x2−I0w¨0δw0∂Mxx∂x−Qxz−I1u¨0−I2θ¨δθ∫−h/2h/2∂Dx∂xcosαz+DzαsinαzδΦdx+Nxxδu0+Mxxδθ+Qxz+NxxT+NxxE∂w0∂xδw0+∫−h/2h/2DxcosαzδΦdz|0L.According to the first term of equation (19) and the basic principle of variations, the coefficients of δu0, δw0, δθ, and δΦ should be set to zero and can be written as(20)δu0:∂Nxx∂x=I0u¨0+I1θ¨,δw0:∂Qxz∂x−NxxT+NxxE∂2w0∂x2=I0w¨0,δθ:∂Mxx∂x−Qxz=2I1u¨0+I2θ¨,δΦ:∫−h/2h/2∂Dx∂xcosαz+Dzαsinαzdz=0.The coordinate boundary conditions of the sandwich piezoelectric beam are determined from the second term of (19). The ends of the beam combined with the corresponding spring constraint can be expressed as x=0,L:(21)Nxx=±kxu0,Qxz+NxxT+NxxE∂w0∂x=±kzw0,Mxx=±kθθ,∫−h/2h/2Dxcosαzdz=±kΦΦ.Therefore, setting different spring stiffness values is equivalent to applying different forces at both ends. Clamped (C) and free (F) correspond to a large value ofk and 0, respectively.


### 2.4. Wave Solutions

The Laplace transformation of (19) is carried out to promote the MRRM in Section 4, and the governing equations in the frequency domain are shown as follows:(22)L11L12L13L14L21L22L23L24L31L32L33L34L41L42L43L44u0w0θΦ=0,where the expressions of Liji,j=1,2,3,4 can be calculated with the following:(23)L11=A11∂2∂x2−I0s2,L12=0,L13=B11∂2∂x2−I1s2,L14=E31∂∂x,L21=0,L22=A55∂2∂x2−I0s2−Ne+NxxT∂2∂x2,L23=A55∂∂x,L24=−E15∂2∂x2,L31=B11∂2∂x2−I1s2,L32=−A55∂∂x,L33=D11∂2∂x2−I2s2−A55,L34=F31+E15∂∂x,L41=E31∂∂x,L42=E15∂2∂x2,L43=F31+E15∂∂x,L44=R11∂2∂x2−R33.In general, the displacement variables that satisfy the governing equations can be set as(24)u0w0θΦ=UeλxWeλxΘeλxΛeλx,where λ is the characteristic wave number and U, W, Θ, and Λ are the wave contribution factors. Submitting equation (24) into equation (22), a set of homogeneous functions about wave amplitude is obtained as follows:(25)T11T12T13T14T21T22T23T24T31T32T33T34T41T42T43T44UWΘΛ=0,where the expressions of Tij can be written as follows:(26)T11=A11λ2−I0s2,T12=0,T13=B11λ2−I1s2,T14=E31λ,T21=0,T22=A55λ2−I0s2−Ne+NxxTλ2,T23=A55λ,T24=−E15λ2,T31=B11λ2−I1s2,T32=−A55λ,T33=D11λ2−I2s2−A55,T34=F31+E15λ,T41=E31λ,T42=E15λ,T43=F31+E15λ,T44=R11λ2−R33.For the general solution of (25), we can obtain the characteristic wave number λ by setting T as zero:(27)b8λ8+b6λ6+b4λ4+b2λ2+b0=0,where b8,b6,b4,b2,b0 are the coefficients corresponding to the wave numbers determined by matrix T. The four pairs of wave number ±λii=1,2,3,4 and corresponding eight eigenvectors can be obtained by equation (27). According to the general method of solving matrix equations, the elements in UWΘΛT can be solved by making one element equal to 1. Thus, the basic solution for the displacement variables can be defined as follows:(28)αi,βi,γi,1,where the coefficients αi,βi,andγi are given by the following:(29)αi=−T14T12T13−T24T22T23−T34T32T33T11T12T13T21T22T23T31T32T33|λ=λi,βi=T11−T14T13T21−T24T23T31−T34T33T11T12T13T21T22T23T31T32T33|λ=λi,γi=T11T12−T14T21T22−T24T31T32−T34T11T12T13T21T22T23T31T32T33|λ=λi.Based on equation (28) and the superposition principle, the displacement variables can be written as follows:(30)u0x=∑i=14αaiaieλix+αdidie−λix,w0x=∑i=14βaiaieλix+βdidie−λix,θx=∑i=14γaiaieλix+γdidie−λix,Φx=∑i=14aieλix+die−λix,where ai=i=1,2,3,4 are the amplitude of the arriving waves and di=i=1,2,3,4 are the departing waves.


## 2.1. Model Description

Figure1 shows a sandwich piezoelectric beam with piezoelectric sheets as the top and bottom layers and a general metal material layer in the middle. The length and thickness are set as L and h, respectively. The piezoelectric and metal layers have a thickness of hp and hm, respectively. Each layer is assumed to be perfectly pasted, and the influence of the adhesive layer is disregarded. For simulation of the boundary force at both ends of the beam, there are three sets of mechanical springs, including two sets of displacement springs and a set of rotary springs, which are set in each group. Furthermore, a set of electric displacement constraints is considered. In the following derivation, the piezoelectric material equation will be described first. For general metal layers, only the piezoelectric term in the equation needs to be removed. Moreover, zk + 1 and zk indicate the distance from the top and bottom of the kth layer to the middle face of the beam, respectively.Figure 1
Piezoelectric composite laminated beam with elastic boundary.
(a)(b)


## 2.2. Kinematic Equation

The external electric potentialϕ and temperature rise ΔT are introduced to describe the thermo-electro-elastic coupling equation of the sandwich piezoelectric beam. According to the assumption of the FSDT, the displacement of any point can be written as follows:(1)uxx,z,t=u0x,t+zθx,t,uzx,z,t=w0x,t.In (1), u0,w0,θ are the displacements and rotation angle of a node in the middle plane. The variables of external electric potential are given by the following [32]:(2)Φ¯x,z,t=−cosαzΦx,t+2zΦ0h,where α=π/h and Φ0 is the initial external electric potential.According to (1) and (2), the mechanical strains and electric field are, respectively, obtained as follows:(3)εxx=∂u∂x,=∂u0∂x+z∂θ∂x,γxz=∂w∂x+∂u∂z,=∂w0∂x+θ,Ex=−∂Φ¯∂x,=cosαz∂Φ∂x,Ez=−∂Φ¯∂z,=−αsinαzΦ−2Φ0h.Based on the stress-strain relationship of orthotropic materials and thermal stress, the constitutive equations of thekth layer from the sandwich piezoelectric beam can be written as(4)σxxk=c¯11kεxx−e¯31kEz−c¯11kα1kΔT,τxzk=c¯55kγxz−e¯15kEx,Dxk=e¯15kγxz+s¯11kEx,Dzk=e¯31kεxx+s¯33kEz+p¯3kΔT,where the overline indicates the transformed coefficient. c¯ij is the elastic constants; e¯ij and s¯ij are piezoelectric and dielectric constants, respectively. Notably, e¯ij and s¯ij are equal to zero when the kth layer is a general metal material; these constants can be specifically expressed as follows [23]:(5)c¯11=c11−c132c33,c¯55=c55,s¯11=s11,s¯33=s33+e332c33,e¯15=e15,e¯31=e31−e33c13c33,p¯3=p3−e33p3c33,where σxx and τxy represent the normal and shear stresses, respectively. Dx and Dz represent the electric displacement. Dy and Ey are neglected, and σyy=σzz=τxy=τyz=0. α1 and p3 correspond to the thermal expansion and pyroelectric constants, respectively. Thus, the thermal strain can be written as(6)εxxTk=α¯1kΔT.By integrating the mechanical stress in (4) in the thickness direction, the resultant expression of stretching stress, moment, and transverse force can be calculated with the following:(7)Nxx=∑k=1n∫zkzk+1σxxkdz,Mxx=∑k=1n∫zkzk+1zσxxkdz,Qxz=∑k=1n∫zkzk+1τxzkdz.Similarly, the corresponding stiffness coefficients can be obtained in the same integration way:(8)A11,B11,D11=∑k=1n∫zkzk+1c¯11k1,z,z2dz,A55=∑k=1n∫zkzk+1c¯55kdz,E31,F31=∑k=1n∫zkzk+1e¯31k1,zαsinαzdz,E15=∑k=1n∫zkzk+1e¯15kcosαzdz,R11=∑k=1n∫zkzk+1s¯11kcos2αzdzc,R33=∑k=1n∫zkzk+1s¯33kα2sin2αzdz,where κ represents the shear correction factor that is determined by a cross-sectional shape. Aij, Bij, and Dij refer to the mechanical stiffnesses, and Eij and Fij denote the electromechanical stiffnesses. Rij stands for the dielectric stiffness [30]. Therefore, we can obtain the constitutive equations of the piezoelectric laminated beam as follows:(9)Nxx=A11∂u0∂x+B11∂θ∂x+E31Φ,Mxx=B11∂u0∂x+D11∂θ∂x+F31Φ,Qxz=κA55∂w0∂x+θ−κE15∂Φ∂x,∫−h/2h/2Dxcosαz=E15∂w0∂x+θ+R11∂Φ∂x,∫−h/2h/2Dzαsinαz=E31∂u0∂x+F31∂θ∂x−R33Φ.


## 2.3. Governing Equations and Geometric Boundary Conditions

The coupling equation under thermal and electric fields of the sandwich piezoelectric beam is obtained by Hamilton’s principle defined by the following:(10)∫t1t2∫VδUM+UT+UE−KdVdt=0,where δ represents the variation and UM indicates the virtual strain energy and virtual electric internal energy calculated by the following [30]:(11)UM=∑k=1n∫0L∫zkzk+1σxxMεxxM+τxzγxz−DxEx−DzEzdzdx.Substituting equations (3), (4), and (7) into (11) and performing variational operations on the above equation, the expression of total virtual internal work is shown as(12)δUM=∫0LNxx∂δu0∂x+Mxx∂δθ∂x+Qxz∂δw0∂x+δθ∫−h/2h/2−Dxcosαz∂δΦ∂x−DzαsinαzδΦdzdx.UT and UE represent the external thermal potential energy and electric potential energy, respectively. In this article, only the axial strain component should be considered, and it can be expressed as(13)εxxT=12∂w∂x2=12∂w0∂x2.Thus, the external thermal and electric potential energies can be written with the following equations:(14)δUT=−∑k=1nm∫0L∫zkzk+1σxxTkδεxxTkdzdx=−∫0LNxxT∂w0∂xδ∂w0∂xdx,δUE=−∑k=1np∫0L∫zkzk+1DzkδεxxTkdzdx=−∫0LNxxE∂w0∂xδ∂w0∂xdx,where NxxT and NxxE are the normal forces as(15)NxxT=c11kεxxTk,NxxE=−2e¯31ϕ0.From the beam displacement field in (1), K represents the kinetic energy of the sandwich piezoelectric beam:(16)∫t1t2δK=−δ∫t1t2∫0L∫−h/2h/2ρku˙2+w˙2dzdx,where the material density of the kth layer is expressed as pk, and the superscript indicates the integral over time.(17)δK=−∑k=1n∫0L∫zkzk+1ρku¨δu+w¨δwdzdx=∑k=1n∫0Lρku˙02+zθ˙2+2zu˙0θ˙+w˙02dx=∫0LI0u¨0δu0+I0w¨0δw0+I1θ¨δu0+I1u¨0δθ+I2θ¨δθdx,where I0,I1,I2 are the kth layer’s inertia coefficients defined by(18)I0,I1,I2=∑k=1n∫zkzk+1ρk1,z,z2dz.Equations (12), (14), and (17) are substituted into equation (10) to derive the governing equations of the sandwich piezoelectric beam under the coupling environment as follows:(19)0=−∫0L∂Nxx∂x−I0u¨0−I1θ¨δu0∂Qxz∂x−NxxT+NxxE∂2w0∂x2−I0w¨0δw0∂Mxx∂x−Qxz−I1u¨0−I2θ¨δθ∫−h/2h/2∂Dx∂xcosαz+DzαsinαzδΦdx+Nxxδu0+Mxxδθ+Qxz+NxxT+NxxE∂w0∂xδw0+∫−h/2h/2DxcosαzδΦdz|0L.According to the first term of equation (19) and the basic principle of variations, the coefficients of δu0, δw0, δθ, and δΦ should be set to zero and can be written as(20)δu0:∂Nxx∂x=I0u¨0+I1θ¨,δw0:∂Qxz∂x−NxxT+NxxE∂2w0∂x2=I0w¨0,δθ:∂Mxx∂x−Qxz=2I1u¨0+I2θ¨,δΦ:∫−h/2h/2∂Dx∂xcosαz+Dzαsinαzdz=0.The coordinate boundary conditions of the sandwich piezoelectric beam are determined from the second term of (19). The ends of the beam combined with the corresponding spring constraint can be expressed as x=0,L:(21)Nxx=±kxu0,Qxz+NxxT+NxxE∂w0∂x=±kzw0,Mxx=±kθθ,∫−h/2h/2Dxcosαzdz=±kΦΦ.Therefore, setting different spring stiffness values is equivalent to applying different forces at both ends. Clamped (C) and free (F) correspond to a large value ofk and 0, respectively.


## 2.4. Wave Solutions

The Laplace transformation of (19) is carried out to promote the MRRM in Section 4, and the governing equations in the frequency domain are shown as follows:(22)L11L12L13L14L21L22L23L24L31L32L33L34L41L42L43L44u0w0θΦ=0,where the expressions of Liji,j=1,2,3,4 can be calculated with the following:(23)L11=A11∂2∂x2−I0s2,L12=0,L13=B11∂2∂x2−I1s2,L14=E31∂∂x,L21=0,L22=A55∂2∂x2−I0s2−Ne+NxxT∂2∂x2,L23=A55∂∂x,L24=−E15∂2∂x2,L31=B11∂2∂x2−I1s2,L32=−A55∂∂x,L33=D11∂2∂x2−I2s2−A55,L34=F31+E15∂∂x,L41=E31∂∂x,L42=E15∂2∂x2,L43=F31+E15∂∂x,L44=R11∂2∂x2−R33.In general, the displacement variables that satisfy the governing equations can be set as(24)u0w0θΦ=UeλxWeλxΘeλxΛeλx,where λ is the characteristic wave number and U, W, Θ, and Λ are the wave contribution factors. Submitting equation (24) into equation (22), a set of homogeneous functions about wave amplitude is obtained as follows:(25)T11T12T13T14T21T22T23T24T31T32T33T34T41T42T43T44UWΘΛ=0,where the expressions of Tij can be written as follows:(26)T11=A11λ2−I0s2,T12=0,T13=B11λ2−I1s2,T14=E31λ,T21=0,T22=A55λ2−I0s2−Ne+NxxTλ2,T23=A55λ,T24=−E15λ2,T31=B11λ2−I1s2,T32=−A55λ,T33=D11λ2−I2s2−A55,T34=F31+E15λ,T41=E31λ,T42=E15λ,T43=F31+E15λ,T44=R11λ2−R33.For the general solution of (25), we can obtain the characteristic wave number λ by setting T as zero:(27)b8λ8+b6λ6+b4λ4+b2λ2+b0=0,where b8,b6,b4,b2,b0 are the coefficients corresponding to the wave numbers determined by matrix T. The four pairs of wave number ±λii=1,2,3,4 and corresponding eight eigenvectors can be obtained by equation (27). According to the general method of solving matrix equations, the elements in UWΘΛT can be solved by making one element equal to 1. Thus, the basic solution for the displacement variables can be defined as follows:(28)αi,βi,γi,1,where the coefficients αi,βi,andγi are given by the following:(29)αi=−T14T12T13−T24T22T23−T34T32T33T11T12T13T21T22T23T31T32T33|λ=λi,βi=T11−T14T13T21−T24T23T31−T34T33T11T12T13T21T22T23T31T32T33|λ=λi,γi=T11T12−T14T21T22−T24T31T32−T34T11T12T13T21T22T23T31T32T33|λ=λi.Based on equation (28) and the superposition principle, the displacement variables can be written as follows:(30)u0x=∑i=14αaiaieλix+αdidie−λix,w0x=∑i=14βaiaieλix+βdidie−λix,θx=∑i=14γaiaieλix+γdidie−λix,Φx=∑i=14aieλix+die−λix,where ai=i=1,2,3,4 are the amplitude of the arriving waves and di=i=1,2,3,4 are the departing waves.


## 3. Solving Method

Before using the MRRM to formulate the structure, the displacement vectorδ and force vector f are introduced first, and equation (30) is rewritten in the following form:(31)δ=u0,w0,θ,ΦT=AδP−xa+DδPxd,f=Nxx,Qxz,Mxx,DxcosαzT=AfP−xa+DfPxd,where(32)a=a1,a2,a3,a4T,d=d1,d2,d3,d4T,Px=diage−λ1x,e−λ2x,e−λ3x,e−λ4xrepresent the arriving wave vectors, departing wave vectors, and phase propagation matrix. Aδ,Dδ and Af,Df are the generalized displacement coefficients and generalized force coefficient matrix, respectively:(33)Aδ=αa1αa2αa3αa4βa1βa2βa3βa4γa1γa2γa3γa41111,Dδ=αd1αd2αd3αd4βd1βd2βd3βd4γd1γd2γd3γd41111,Af=Af11Af12Af13Af14Af21Af22Af23Af24Af31Af32Af33Af34Af41Af42Af43Af44,Df=Df11Df12Df13Df14Df21Df22Df23Df24Df31Df32Df33Df34Df41Df42Df43Df44,where(34)Af1jj=1,2,3,4=A11αajλj−E31+B11γajλj,Af2jj=1,2,3,4=A55γaj−E15λj+A55βajλj−NxxEβajλj2−NxxTβajλj2,Af3jj=1,2,3,4=B11αajλj−F31+D11γajλj,Af4jj=1,2,3,4=E15γaj+R11λj+E15βajλj,Df1jj=1,2,3,4=−A11αajλj−E31−B11γajλj,Df2jj=1,2,3,4=A55γaj+E15λj−A55βajλj−NxxEβajλj2−NxxTβajλj2,Df3jj=1,2,3,4=−B11αajλj−F31−D11γajλj,Df4jj=1,2,3,4=E15γaj−R11λj−E15βajλj.


### 3.1. Generalized Phase Relationship

As the first step in constructing the MRRM formulation, a pair of local dual coordinate systems should be established for each unit on the structure (Figure2). Nodes 1 and 2 are considered the two ends of the structure, and the positive direction is set from node 1 to node 2 as x12. The coordinate systems at nodes 1 and 2 are defined as O1−X12Z12 and O2−X21Z21, respectively. If the first coordinate system O1−X12Z12 is defined to be parallel with the positive direction, then the stiffness coefficient and inertia moment of the structure in the “21” coordinate system will change accordingly. The related coefficients in the “21” coordinate system can be written as follows:(35)A1121=A1112,B1121=−B1112,D1121=D1112,E3121=E3112,E1521=E1512,F3121=−F3112,R1121=R1112,R3321=R3312,NxxT,21=NxxT,12,I021=I012,I121=−I112,I221=I212.Figure 2
Dual coordinate system of the composite laminated beam.Therefore, the relations of the generalized displacement coefficientsAδ,Dδ and Af,Df in the dual coordinate system can be obtained as follows:(36)Aδ21=−TδDδ12,Dδ21=−TδAδ12,Af21=−TfDf12,Df21=−TfAf12,where matrix Tδ=diag−1,−1,1,1 and Tf=diag1,1,−1,1. Although two wave solutions can be obtained, the physical quantity at the same point on the beam is constant, and the vectors δ and f in the dual coordinate system are expressed with the following:(37)δ12x12=Tδδ21L−x12,f12x12=Tff21L−x12.Before deriving the phase relationship of the beam structure, the wave solutions need to be classified first. When the wave direction is the same or opposite withx12, it is called the departing or arriving wave, respectively. Thus, substituting equation (37) into equation (31), the relationship between aij and dij can be simplified and classified as follows:(38)Aδ21P0a21+Dδ21P0d21=TδAδ12P−La12+Dδ12PLd12,Af21P0a21+Df21P0d21=TfAf12P−La12+Df12PLd12,(39)a12a21=Phd12d21,where Ph is the phase matrix of the sandwich piezoelectric beam.


### 3.2. Scattering Relation

Using only the phase relationship of the structure is insufficient to obtain the dynamic characteristics of the sandwich piezoelectric beams. The scattering relationship at each node can be obtained through geometric equations and force balance conditions. Given that the spring boundary exists at nodes 1 and 2 on both ends, equation (31) can be used to establish the relationship between generalized displacements and forces through the spring stiffness matrix:(40)K1,2=diagku0,kw0,kθ,kΦ.Thus, we can obtain the following:(41)K1Aδ12a12+Dδ12d12=Af12a12+Df12d12,K2Aδ21a21+Dδ21d21=Af21a21+Df21d21.Table1 provides the stiffness matrix including the C, S, F, and elastic support (E) for common boundary conditions in engineering. The value of 1018 represents a very large stiffness coefficient to achieve the C condition because an infinite value cannot be calculated. The scattering relation can be obtained by classifying vectors aij and dij in the same manner:(42)d12d21=Sa12a21,where S represents the scattering matrix of the sandwich piezoelectric beam.Table 1
Boundary conditions and stiffness matrixes.
Boundary conditionsEssential conditionsStiffness matrixesFNxx = Qxz = Mxx = 0diag (0, 0, 0, 0)Su0 = w0 = Mxx = 0diag (1018, 1018, 0, 0)Cu0 = w0 = θ = 0diag (1018, 1018, 1018, 0)Eu0, w0, θ ≠ 0diag (ku, kw, kθ, kΦ)E1u0 = θ = 0, w0≠0diag (1018, kw, 1018, 0)E2u0 = w0 = 0, θ ≠ 0diag (1018, 1018, kθ, 0)E3u0 = 0, w0, θ ≠ 0diag (1018, kw, kθ, 0)


### 3.3. Reverberation-Ray Matrix and Free Vibration

According to the relation of equations (39) and (42), the overall MRRM equation of the piezoelectric laminated beam without considering external excitation can be obtained as follows:(43)E−Rd=0,where R=SPh and E represents the identity matrix. R is called the reverberation-ray matrix of the present sandwich piezoelectric beam. For free vibration analysis, the determinant of equation (43) is set to zero:(44)detE−Rω=0,where the solution of equation (44) is the natural vibration parameter of the sandwich piezoelectric beam.Equation (44) is a complex transcendental equation about frequency ω and is not easily solved by general analytical methods. However, for this complex equation, as ω gradually increases, the modulus of the complex determinant always assumes a minimum value, which is zero when it passes through its zero point. Therefore, the golden-section search (GSS) method can be used to search for natural frequencies. The evaluation function is defined as follows:(45)fω=detE−Rw.The specific process is performed as follows:Step 1: define the initial trial frequencyω0, frequency search step size Δω, and maximum frequency ωmax∘Step 2: generating two search frequenciesω1=ω0+Δω and ω2=ω1+Δω, calculate the evaluation function values corresponding to the three frequencies f0,f1,andf2Step 3: iff0<f1 and f1<f2, proceed to Step 4; else, return to Step 2 with resetting ω0=ω1Step 4: execute the GSS method to search the undetermined frequencyω3, and calculate the value of function f3Step 5: when the maximum frequencyωmax is reached, stop the search; else, return to Step 2 to reset ω0=ω3+ΔωFigure3 shows the specific flowchart of the GSS algorithm.Figure 3
Natural frequency calculation process based on the golden-section search algorithm.


## 3.1. Generalized Phase Relationship

As the first step in constructing the MRRM formulation, a pair of local dual coordinate systems should be established for each unit on the structure (Figure2). Nodes 1 and 2 are considered the two ends of the structure, and the positive direction is set from node 1 to node 2 as x12. The coordinate systems at nodes 1 and 2 are defined as O1−X12Z12 and O2−X21Z21, respectively. If the first coordinate system O1−X12Z12 is defined to be parallel with the positive direction, then the stiffness coefficient and inertia moment of the structure in the “21” coordinate system will change accordingly. The related coefficients in the “21” coordinate system can be written as follows:(35)A1121=A1112,B1121=−B1112,D1121=D1112,E3121=E3112,E1521=E1512,F3121=−F3112,R1121=R1112,R3321=R3312,NxxT,21=NxxT,12,I021=I012,I121=−I112,I221=I212.Figure 2
Dual coordinate system of the composite laminated beam.Therefore, the relations of the generalized displacement coefficientsAδ,Dδ and Af,Df in the dual coordinate system can be obtained as follows:(36)Aδ21=−TδDδ12,Dδ21=−TδAδ12,Af21=−TfDf12,Df21=−TfAf12,where matrix Tδ=diag−1,−1,1,1 and Tf=diag1,1,−1,1. Although two wave solutions can be obtained, the physical quantity at the same point on the beam is constant, and the vectors δ and f in the dual coordinate system are expressed with the following:(37)δ12x12=Tδδ21L−x12,f12x12=Tff21L−x12.Before deriving the phase relationship of the beam structure, the wave solutions need to be classified first. When the wave direction is the same or opposite withx12, it is called the departing or arriving wave, respectively. Thus, substituting equation (37) into equation (31), the relationship between aij and dij can be simplified and classified as follows:(38)Aδ21P0a21+Dδ21P0d21=TδAδ12P−La12+Dδ12PLd12,Af21P0a21+Df21P0d21=TfAf12P−La12+Df12PLd12,(39)a12a21=Phd12d21,where Ph is the phase matrix of the sandwich piezoelectric beam.


## 3.2. Scattering Relation

Using only the phase relationship of the structure is insufficient to obtain the dynamic characteristics of the sandwich piezoelectric beams. The scattering relationship at each node can be obtained through geometric equations and force balance conditions. Given that the spring boundary exists at nodes 1 and 2 on both ends, equation (31) can be used to establish the relationship between generalized displacements and forces through the spring stiffness matrix:(40)K1,2=diagku0,kw0,kθ,kΦ.Thus, we can obtain the following:(41)K1Aδ12a12+Dδ12d12=Af12a12+Df12d12,K2Aδ21a21+Dδ21d21=Af21a21+Df21d21.Table1 provides the stiffness matrix including the C, S, F, and elastic support (E) for common boundary conditions in engineering. The value of 1018 represents a very large stiffness coefficient to achieve the C condition because an infinite value cannot be calculated. The scattering relation can be obtained by classifying vectors aij and dij in the same manner:(42)d12d21=Sa12a21,where S represents the scattering matrix of the sandwich piezoelectric beam.Table 1
Boundary conditions and stiffness matrixes.
Boundary conditionsEssential conditionsStiffness matrixesFNxx = Qxz = Mxx = 0diag (0, 0, 0, 0)Su0 = w0 = Mxx = 0diag (1018, 1018, 0, 0)Cu0 = w0 = θ = 0diag (1018, 1018, 1018, 0)Eu0, w0, θ ≠ 0diag (ku, kw, kθ, kΦ)E1u0 = θ = 0, w0≠0diag (1018, kw, 1018, 0)E2u0 = w0 = 0, θ ≠ 0diag (1018, 1018, kθ, 0)E3u0 = 0, w0, θ ≠ 0diag (1018, kw, kθ, 0)


## 3.3. Reverberation-Ray Matrix and Free Vibration

According to the relation of equations (39) and (42), the overall MRRM equation of the piezoelectric laminated beam without considering external excitation can be obtained as follows:(43)E−Rd=0,where R=SPh and E represents the identity matrix. R is called the reverberation-ray matrix of the present sandwich piezoelectric beam. For free vibration analysis, the determinant of equation (43) is set to zero:(44)detE−Rω=0,where the solution of equation (44) is the natural vibration parameter of the sandwich piezoelectric beam.Equation (44) is a complex transcendental equation about frequency ω and is not easily solved by general analytical methods. However, for this complex equation, as ω gradually increases, the modulus of the complex determinant always assumes a minimum value, which is zero when it passes through its zero point. Therefore, the golden-section search (GSS) method can be used to search for natural frequencies. The evaluation function is defined as follows:(45)fω=detE−Rw.The specific process is performed as follows:Step 1: define the initial trial frequencyω0, frequency search step size Δω, and maximum frequency ωmax∘Step 2: generating two search frequenciesω1=ω0+Δω and ω2=ω1+Δω, calculate the evaluation function values corresponding to the three frequencies f0,f1,andf2Step 3: iff0<f1 and f1<f2, proceed to Step 4; else, return to Step 2 with resetting ω0=ω1Step 4: execute the GSS method to search the undetermined frequencyω3, and calculate the value of function f3Step 5: when the maximum frequencyωmax is reached, stop the search; else, return to Step 2 to reset ω0=ω3+ΔωFigure3 shows the specific flowchart of the GSS algorithm.Figure 3
Natural frequency calculation process based on the golden-section search algorithm.


## 4. Result Verification and Discussion

### 4.1. Result Verification

Currently, several types of numerical examples are used to study the free vibration characteristics of sandwich piezoelectric beams. First, the calculated results are compared with the findings of literature studies and software simulations to verify the calculation correctness. On this basis, a series of parameter coupling studies, including elastic modulus, ply orientation of materials, thickness ratio, external electric potential, and temperature rise, is considered. The effect of single and coupling factors on vibration characteristics will be investigated separately. Table2 lists the relevant material parameters for verification and subsequent parameterization studies. The nondimensional results Ω=ωLρ/c55PZT−5A will be used to express the frequency parameters if no extra instructions are given in the following examples.Table 2
Related material parameters.
PZT-5A [55, 56]c11 = 121 GPa, c13 = 75.2 GPa, c33 = 111 GPa, c55 = 21.1 GPa, ρ = 7750 kg/m3, α1 = 2.68 × 10−6/Ke31 = −5.4 C/m2, e33 = 15.8 C/m2, e15 = 12.3 C/m2s11 = 8.107 nF/m, s33 = 7.34 nF/mAl [56]c11 = 72.2 GPa, c55 = 26.94 GPa, ρ = 2700 kg/m3, α1 = 23.2 × 10−6/KMaterial-Ac11 = 138.4 GPa, c55 = 162.5 GPa, ρ = 7500 kg/m3Material-Bc11 = 157.3 GPa, c55 = 45.5 GPa, ρ = 5676 kg/m3Material-CE = open, μ = 0.3, ρ = 2700 kg/m3Material-DE1/E2 = 15, E2 = E3, μ12 = μ13 = μ23 = 0.25, G12 = G13 = 0.6E2, G23 = 0.5E2, ρ = 1500 kg/m3, α1 = 6 × 10−6/KFirst, a double-layer composite beam composed of materials A and B is introduced without the thermal environment and piezoelectric effect. As mentioned earlier, the general material layer only needs to set the piezoelectric coefficients in the equation to zero. He et al. [32] analyzed this example based on the DQM and considered the boundary conditions and slenderness ratio. Table 3 shows the comparison results of the first six dimensionless frequencies under different boundary conditions and slenderness ratios, where ω is defined as Ω=ωhρ/c55A. The numerical errors are in a relatively small range, and the poorest result error is at 1.7%.Table 3
Comparison of natural frequencies for a double-layer laminated beam under different boundary conditions and slenderness ratios (material-A/material-B).
BCSlenderness ratioMode123457C-FL/h = 15[33]0.00470.02880.0790.10990.15020.2395Present0.00470.02930.08020.11030.15260.2434Error (%)0.00−1.71−1.50−0.36−1.57−1.60L/h = 7[33]0.02130.12530.23580.32340.5760.7063Present0.02160.12720.23630.32850.58530.709Error (%)−1.39−1.49−0.21−1.55−1.59−0.38C-CL/h = 15[33]0.02920.07860.14950.22050.23840.3424Present0.02970.07980.15170.22060.24170.3466Error (%)−1.68−1.50−1.45−0.05−1.37−1.21L/h = 7[33]0.12640.31840.47190.56750.85150.9425Present0.12820.32280.47260.57420.860.9453Error (%)−1.40−1.36−0.15−1.17−0.99−0.30Next, the material in the above numerical example is replaced with piezoelectric material, a double-layer beam with equal thickness composed of PZT-5A analyzed by Fernandes and Pouget [11]; in consideration of the Bernoulli–Euler beam theory, a purely mechanical beam theory and an electromechanical coupling model under an electric field environment are established, and the related material parameters can be found in Table 2 (dimensions: L = 0.4 m and h = 0.01 m). Table 4 lists the first six dimensionless frequencies in different mechanical boundary conditions, which show that our frequency parameters are notably close to those of Michael Krommer and FEM simulation results. In the table, “NE” stands for the layers not equipped with electrodes, which are consistent with the model in this example.Table 4
Comparison of natural frequencies for a double-layer piezoelectric laminated beam under different boundary conditions (PZT-5A/PZT-5A,L = 0.4 m, and h = 0.01 m).
BCMode123456S-SNE [13]0.12960.51851.16652.07393.24034.6662FEM0.12970.51761.16002.05073.18154.5433Present0.12960.51741.15932.04923.17864.5387Error (%)0.080.040.060.070.090.10C-CNE [13]0.29380.80991.58782.62463.92095.4763FEM0.29420.80711.57242.57943.81875.2791Present0.29320.80491.56862.57383.81115.2691Error (%)0.340.270.240.220.200.19Finally, a sandwich piezoelectric laminated beam withL = 0.4 m and h = 0.03 m with equal thickness of each layer composed of PZT-5A/aluminum/PZT-5A is considered. This verification example and the following parameter studies will all use this model unless special instructions are specified. In the comparison with Ansys v18.2 software, 80 solid 5-grid elements and 40 solid 45-grid elements are used to simulate the piezoelectric and metal layers, respectively. The external coefficients Φ0 and △T are set to 10 kV and 25–200°C. Table 5 shows the first six dimensionless natural frequencies compared with Ansys v18.2 software under the clamped-clamped boundary condition. The frequency values of the two variables are extremely close, and the largest error is 0.47%.Table 5
Comparison of natural frequencies for a sandwich piezoelectric laminated beam under different temperature rises (PZT-5A/Al/PZT-5A,L = 0.4 m, h = 0.03 m, C-C, and Φ0 = 10 kV).
T (°C)Mode12345625Present0.96062.55294.77507.478910.549313.8948FEM0.96432.55524.76807.457710.516013.8619Error (%)0.390.09−0.15−0.28−0.31−0.2450Present0.95422.54414.76517.468210.537813.8827FEM0.95792.54644.75817.446710.504213.8491Error (%)0.390.09−0.15−0.29−0.32−0.24100Present0.94102.52634.74517.446710.515013.8584FEM0.94502.52854.73807.424810.480313.8233Error (%)0.430.09−0.15−0.29−0.33−0.25200Present0.91412.49004.70517.403610.469113.8096FEM0.91842.49244.69767.380610.432613.7720Error (%)0.470.10−0.16−0.31−0.35−0.27Based on the above comparisons and verifications, the proposed method is a simple and effective technique for the analysis of the vibration characteristics of piezoelectric composite beams, and it provides a convenient adjustment interface for various boundary conditions and coupling of different external conditions or parameters. In this manner, the following parameterization studies require minor changes to improve the computational efficiency and accuracy.


### 4.2. Case 1: Elastic Modulus and Boundary Conditions in the Thermal Environment

Different from general composite materials or functionally graded beams, the sandwich piezoelectric laminated structure has a good design capability because it is composed of different kinds of materials, and various parameters have varied effects on its dynamic performance. This paper mainly studies five aspects of variables, including temperature rise, external electric potential, material and geometric properties of the middle layer, elastic boundary conditions, and the free vibration characteristic under coupled fields.For the detailed investigation of the effect of spring restraint constant on frequency characteristics, three elastic boundary conditions (E1, E2, and E3) are additionally defined, and they correspond to thez-direction spring constant, rotational spring constant, and condition when both constants are considered. Table 1 provides the detailed essential conditions and spring coefficients. With the calculated natural frequency under different boundary combinations, the changing rules in Figure 4 show that the frequency parameters of each order increase within a certain elastic coefficient interval that is generally between 105 and 1010. The interval gradually moves to the right as the frequency order increases. Furthermore, different boundary combinations will result in different frequency characteristics. Figures 4(a) and 4(b) show the independent influence of kw and kθ when one end is C or S, respectively. Figure 4(d) shows the effect of kw and kθ acting at the same time when one end is C or S. Figures 4(c) and 4(e) show that the effect of both ends is elastic restraint.Figure 4
Variations of the frequency parameter Ω versus the elastic restraint parameters for a sandwich laminated beam (PZT-5A/Al/PZT-5A).
(a)(b)(c)(d)(e)Next, the effects of elastic boundary and traditional boundary conditions on frequency parameters are studied by giving a specific spring constant. Based on the spring constants corresponding to the frequency variation interval in Figure4, two defined elastic boundary conditions Ea and Eb are considered besides the conventional S, C, and F conditions, where Ea indicates that only the z-direction is elastically restrained (w0≠0), and the spring constant kw is set as 106; Eb indicates that the rotation direction and z-direction are elastically restrained (w0, θ≠0), and kw and kθ are set as 106. Figures 5(a)–5(d) show the frequency parameters under different boundary conditions. Figures 5(c) and 5(d) illustrate that the frequency parameters can be further adjusted by changing the spring constant.Figure 5
Variations of the frequency parameter Ω versus boundary conditions and elastic modulus under different temperature rises for a sandwich laminated beam (PZT-5A/material-C/PZT-5A, material-C,E = open, α1 = 23.2 × 10−6/K). (a) C-C. (b) S-S. (c) Ea-Eb. (d) Ea-Ea.
(a)(b)(c)(d)On the basis of the above studies, we further consider the effect of the properties of the middle metal layer under the thermal environment on the frequency parameters. First, the elastic modulus of the middle layer is set from 30 GPa to 350 GPa. However, the density and Poisson’s ratio remain unchanged, and the range of temperature rise is 0°C–200°C. Figure5 also shows the change trend of frequency parameters under different elastic modulus and temperature rise conditions. The thermal environment has a negative effect on the free vibration of the piezoelectric beam because the stiffness is weakened when the temperature rise causes internal stress in the beam. From the horizontal view in Figure 5, when Poisson’s ratio and density of the core layer are constant, the free vibration parameters will increase with the increase of its elastic modulus, and high-order frequencies are more evident than the low-order ones.


### 4.3. Case 2: Ply Orientation and the Number of Layers under Different Electric Potentials

In this section, the middle metal layer is regarded as a composite layer composed of anisotropic materials (Figure6(a)). The effects of the lamination method, ply orientation, the number of layers, and external electric potential on the frequency parameters are discussed. Table 2 gives the parameters of material-D, and there are five types of lamination method, denoted as θ, 0/θ, θ/0/θ, 0/θ/0/θ, and θ/0/θ/0/θ, where θ represents the ply orientation to the x-direction which is set as 0°–180°. The external electric potential is set from −20 kV to 20 kV. Figure 7 shows the frequency parameters of the clamped-clamped (C-C) beam under the five types of lamination methods. Figure 7(a) shows that the frequency parameters are symmetrical at about x = 90°. For different lamination methods, the frequency range of beams with more layers is significantly larger than those with few layers. However, the overall trend is similar. From the horizontal view, a ply orientation with the x-axis has a negative effect on free vibration parameters, and the slope is large between 30°–60° and 120°–150°. Figure 7(b) shows the effect of ply orientation under different lamination methods. When the size of the laminated beam is constant and the angle is zero, the effect of the number of middle layers on free vibration can be obtained. It is obvious that more middle layers indicate high-frequency parameters. When the ply orientation becomes larger, the change mode of different lamination methods varies although the overall trend of the frequency is still increasing. The comparison reveals that it is related to the proportion of the zero-degree layers to the total layers, and the range of frequency changes decreases as this ratio increases. In Figure 7(c), external electric potential is similar with temperature rise that has a negative effect on free vibration parameters. However, the overall effect is notably small.Figure 6
Thickness distribution under different lamination methods.
(a)(b)(c)Figure 7
Variations of the frequency parameter Ω versus the ply orientation, lamination method, and electric potential for a sandwich laminated beam (PZT-5A/material-D/…/PZT-5A, material-D,E1/E2 = 15, E2 = 109 GPa).
(a)(b)(c)


### 4.4. Case 3: Thickness Ratio and Composition Method in the Thermal Environment

Next, we study the frequency characteristics of the thickness ratio under different material combinations in two cases. Similar to Case 1, the middle layer is composed of isotropic materials, and the beam size is constant. As shown in Figure3(b), six thickness ratios (hm/hp = 0.5, 1, 2, 5, 20, and 50) of the metal layer to the piezoelectric layer are considered by calculating the frequency parameters under different thickness ratios and the number of layers under C-C boundary conditions. Figure 8 illustrates that a piezoelectric laminated beam with a large thickness ratio has a high frequency at any number of layers. However, the growth rate decreases as the thickness ratio increases and becomes more evident with the increase in middle layers. Furthermore, if the middle layer is composed of isotropic and anisotropic layers with hm and hc for each layer as shown in Figure 3(c), then the effect of different numbers of layers and thickness ratios (hm/hc) on the frequency parameters should be reconsidered, and the temperature rise is set from 0°C to 200°C. Figure 9(a) shows that when piezoelectric, anisotropic material, and isotropic material layers together form a symmetrical composite structure, the frequency characteristics differ from those of previous examples. For each order frequency, the thickness ratio has a positive effect on free vibration parameters, and the value is between 1 and 5. If the thickness ratio is less than 1 and greater than 5, the frequency will change in the opposite direction. In addition, the change trends in Figures 9(c) and 9(d) are similar with the rules in Case 1. The temperature rise still has a negative but negligible effect on the frequency parameters in this example.Figure 8
Variations of the frequency parameter Ω versus the thickness ratio (hm/hp) and the number of layers for an N-layer piezoelectric laminated beam (PZT-5A/Al/Al/…/PZT-5A).
(a)(b)Figure 9
Variations of the frequency parameter Ω versus the thickness ratio (hm/hc), number of layers, and temperature rise for an N-layer piezoelectric laminated beam (PZT-5A/material-D/Al/…/PZT-5A).
(a)(b)(c)(d)


## 4.1. Result Verification

Currently, several types of numerical examples are used to study the free vibration characteristics of sandwich piezoelectric beams. First, the calculated results are compared with the findings of literature studies and software simulations to verify the calculation correctness. On this basis, a series of parameter coupling studies, including elastic modulus, ply orientation of materials, thickness ratio, external electric potential, and temperature rise, is considered. The effect of single and coupling factors on vibration characteristics will be investigated separately. Table2 lists the relevant material parameters for verification and subsequent parameterization studies. The nondimensional results Ω=ωLρ/c55PZT−5A will be used to express the frequency parameters if no extra instructions are given in the following examples.Table 2
Related material parameters.
PZT-5A [55, 56]c11 = 121 GPa, c13 = 75.2 GPa, c33 = 111 GPa, c55 = 21.1 GPa, ρ = 7750 kg/m3, α1 = 2.68 × 10−6/Ke31 = −5.4 C/m2, e33 = 15.8 C/m2, e15 = 12.3 C/m2s11 = 8.107 nF/m, s33 = 7.34 nF/mAl [56]c11 = 72.2 GPa, c55 = 26.94 GPa, ρ = 2700 kg/m3, α1 = 23.2 × 10−6/KMaterial-Ac11 = 138.4 GPa, c55 = 162.5 GPa, ρ = 7500 kg/m3Material-Bc11 = 157.3 GPa, c55 = 45.5 GPa, ρ = 5676 kg/m3Material-CE = open, μ = 0.3, ρ = 2700 kg/m3Material-DE1/E2 = 15, E2 = E3, μ12 = μ13 = μ23 = 0.25, G12 = G13 = 0.6E2, G23 = 0.5E2, ρ = 1500 kg/m3, α1 = 6 × 10−6/KFirst, a double-layer composite beam composed of materials A and B is introduced without the thermal environment and piezoelectric effect. As mentioned earlier, the general material layer only needs to set the piezoelectric coefficients in the equation to zero. He et al. [32] analyzed this example based on the DQM and considered the boundary conditions and slenderness ratio. Table 3 shows the comparison results of the first six dimensionless frequencies under different boundary conditions and slenderness ratios, where ω is defined as Ω=ωhρ/c55A. The numerical errors are in a relatively small range, and the poorest result error is at 1.7%.Table 3
Comparison of natural frequencies for a double-layer laminated beam under different boundary conditions and slenderness ratios (material-A/material-B).
BCSlenderness ratioMode123457C-FL/h = 15[33]0.00470.02880.0790.10990.15020.2395Present0.00470.02930.08020.11030.15260.2434Error (%)0.00−1.71−1.50−0.36−1.57−1.60L/h = 7[33]0.02130.12530.23580.32340.5760.7063Present0.02160.12720.23630.32850.58530.709Error (%)−1.39−1.49−0.21−1.55−1.59−0.38C-CL/h = 15[33]0.02920.07860.14950.22050.23840.3424Present0.02970.07980.15170.22060.24170.3466Error (%)−1.68−1.50−1.45−0.05−1.37−1.21L/h = 7[33]0.12640.31840.47190.56750.85150.9425Present0.12820.32280.47260.57420.860.9453Error (%)−1.40−1.36−0.15−1.17−0.99−0.30Next, the material in the above numerical example is replaced with piezoelectric material, a double-layer beam with equal thickness composed of PZT-5A analyzed by Fernandes and Pouget [11]; in consideration of the Bernoulli–Euler beam theory, a purely mechanical beam theory and an electromechanical coupling model under an electric field environment are established, and the related material parameters can be found in Table 2 (dimensions: L = 0.4 m and h = 0.01 m). Table 4 lists the first six dimensionless frequencies in different mechanical boundary conditions, which show that our frequency parameters are notably close to those of Michael Krommer and FEM simulation results. In the table, “NE” stands for the layers not equipped with electrodes, which are consistent with the model in this example.Table 4
Comparison of natural frequencies for a double-layer piezoelectric laminated beam under different boundary conditions (PZT-5A/PZT-5A,L = 0.4 m, and h = 0.01 m).
BCMode123456S-SNE [13]0.12960.51851.16652.07393.24034.6662FEM0.12970.51761.16002.05073.18154.5433Present0.12960.51741.15932.04923.17864.5387Error (%)0.080.040.060.070.090.10C-CNE [13]0.29380.80991.58782.62463.92095.4763FEM0.29420.80711.57242.57943.81875.2791Present0.29320.80491.56862.57383.81115.2691Error (%)0.340.270.240.220.200.19Finally, a sandwich piezoelectric laminated beam withL = 0.4 m and h = 0.03 m with equal thickness of each layer composed of PZT-5A/aluminum/PZT-5A is considered. This verification example and the following parameter studies will all use this model unless special instructions are specified. In the comparison with Ansys v18.2 software, 80 solid 5-grid elements and 40 solid 45-grid elements are used to simulate the piezoelectric and metal layers, respectively. The external coefficients Φ0 and △T are set to 10 kV and 25–200°C. Table 5 shows the first six dimensionless natural frequencies compared with Ansys v18.2 software under the clamped-clamped boundary condition. The frequency values of the two variables are extremely close, and the largest error is 0.47%.Table 5
Comparison of natural frequencies for a sandwich piezoelectric laminated beam under different temperature rises (PZT-5A/Al/PZT-5A,L = 0.4 m, h = 0.03 m, C-C, and Φ0 = 10 kV).
T (°C)Mode12345625Present0.96062.55294.77507.478910.549313.8948FEM0.96432.55524.76807.457710.516013.8619Error (%)0.390.09−0.15−0.28−0.31−0.2450Present0.95422.54414.76517.468210.537813.8827FEM0.95792.54644.75817.446710.504213.8491Error (%)0.390.09−0.15−0.29−0.32−0.24100Present0.94102.52634.74517.446710.515013.8584FEM0.94502.52854.73807.424810.480313.8233Error (%)0.430.09−0.15−0.29−0.33−0.25200Present0.91412.49004.70517.403610.469113.8096FEM0.91842.49244.69767.380610.432613.7720Error (%)0.470.10−0.16−0.31−0.35−0.27Based on the above comparisons and verifications, the proposed method is a simple and effective technique for the analysis of the vibration characteristics of piezoelectric composite beams, and it provides a convenient adjustment interface for various boundary conditions and coupling of different external conditions or parameters. In this manner, the following parameterization studies require minor changes to improve the computational efficiency and accuracy.


## 4.2. Case 1: Elastic Modulus and Boundary Conditions in the Thermal Environment

Different from general composite materials or functionally graded beams, the sandwich piezoelectric laminated structure has a good design capability because it is composed of different kinds of materials, and various parameters have varied effects on its dynamic performance. This paper mainly studies five aspects of variables, including temperature rise, external electric potential, material and geometric properties of the middle layer, elastic boundary conditions, and the free vibration characteristic under coupled fields.For the detailed investigation of the effect of spring restraint constant on frequency characteristics, three elastic boundary conditions (E1, E2, and E3) are additionally defined, and they correspond to thez-direction spring constant, rotational spring constant, and condition when both constants are considered. Table 1 provides the detailed essential conditions and spring coefficients. With the calculated natural frequency under different boundary combinations, the changing rules in Figure 4 show that the frequency parameters of each order increase within a certain elastic coefficient interval that is generally between 105 and 1010. The interval gradually moves to the right as the frequency order increases. Furthermore, different boundary combinations will result in different frequency characteristics. Figures 4(a) and 4(b) show the independent influence of kw and kθ when one end is C or S, respectively. Figure 4(d) shows the effect of kw and kθ acting at the same time when one end is C or S. Figures 4(c) and 4(e) show that the effect of both ends is elastic restraint.Figure 4
Variations of the frequency parameter Ω versus the elastic restraint parameters for a sandwich laminated beam (PZT-5A/Al/PZT-5A).
(a)(b)(c)(d)(e)Next, the effects of elastic boundary and traditional boundary conditions on frequency parameters are studied by giving a specific spring constant. Based on the spring constants corresponding to the frequency variation interval in Figure4, two defined elastic boundary conditions Ea and Eb are considered besides the conventional S, C, and F conditions, where Ea indicates that only the z-direction is elastically restrained (w0≠0), and the spring constant kw is set as 106; Eb indicates that the rotation direction and z-direction are elastically restrained (w0, θ≠0), and kw and kθ are set as 106. Figures 5(a)–5(d) show the frequency parameters under different boundary conditions. Figures 5(c) and 5(d) illustrate that the frequency parameters can be further adjusted by changing the spring constant.Figure 5
Variations of the frequency parameter Ω versus boundary conditions and elastic modulus under different temperature rises for a sandwich laminated beam (PZT-5A/material-C/PZT-5A, material-C,E = open, α1 = 23.2 × 10−6/K). (a) C-C. (b) S-S. (c) Ea-Eb. (d) Ea-Ea.
(a)(b)(c)(d)On the basis of the above studies, we further consider the effect of the properties of the middle metal layer under the thermal environment on the frequency parameters. First, the elastic modulus of the middle layer is set from 30 GPa to 350 GPa. However, the density and Poisson’s ratio remain unchanged, and the range of temperature rise is 0°C–200°C. Figure5 also shows the change trend of frequency parameters under different elastic modulus and temperature rise conditions. The thermal environment has a negative effect on the free vibration of the piezoelectric beam because the stiffness is weakened when the temperature rise causes internal stress in the beam. From the horizontal view in Figure 5, when Poisson’s ratio and density of the core layer are constant, the free vibration parameters will increase with the increase of its elastic modulus, and high-order frequencies are more evident than the low-order ones.


## 4.3. Case 2: Ply Orientation and the Number of Layers under Different Electric Potentials

In this section, the middle metal layer is regarded as a composite layer composed of anisotropic materials (Figure6(a)). The effects of the lamination method, ply orientation, the number of layers, and external electric potential on the frequency parameters are discussed. Table 2 gives the parameters of material-D, and there are five types of lamination method, denoted as θ, 0/θ, θ/0/θ, 0/θ/0/θ, and θ/0/θ/0/θ, where θ represents the ply orientation to the x-direction which is set as 0°–180°. The external electric potential is set from −20 kV to 20 kV. Figure 7 shows the frequency parameters of the clamped-clamped (C-C) beam under the five types of lamination methods. Figure 7(a) shows that the frequency parameters are symmetrical at about x = 90°. For different lamination methods, the frequency range of beams with more layers is significantly larger than those with few layers. However, the overall trend is similar. From the horizontal view, a ply orientation with the x-axis has a negative effect on free vibration parameters, and the slope is large between 30°–60° and 120°–150°. Figure 7(b) shows the effect of ply orientation under different lamination methods. When the size of the laminated beam is constant and the angle is zero, the effect of the number of middle layers on free vibration can be obtained. It is obvious that more middle layers indicate high-frequency parameters. When the ply orientation becomes larger, the change mode of different lamination methods varies although the overall trend of the frequency is still increasing. The comparison reveals that it is related to the proportion of the zero-degree layers to the total layers, and the range of frequency changes decreases as this ratio increases. In Figure 7(c), external electric potential is similar with temperature rise that has a negative effect on free vibration parameters. However, the overall effect is notably small.Figure 6
Thickness distribution under different lamination methods.
(a)(b)(c)Figure 7
Variations of the frequency parameter Ω versus the ply orientation, lamination method, and electric potential for a sandwich laminated beam (PZT-5A/material-D/…/PZT-5A, material-D,E1/E2 = 15, E2 = 109 GPa).
(a)(b)(c)


## 4.4. Case 3: Thickness Ratio and Composition Method in the Thermal Environment

Next, we study the frequency characteristics of the thickness ratio under different material combinations in two cases. Similar to Case 1, the middle layer is composed of isotropic materials, and the beam size is constant. As shown in Figure3(b), six thickness ratios (hm/hp = 0.5, 1, 2, 5, 20, and 50) of the metal layer to the piezoelectric layer are considered by calculating the frequency parameters under different thickness ratios and the number of layers under C-C boundary conditions. Figure 8 illustrates that a piezoelectric laminated beam with a large thickness ratio has a high frequency at any number of layers. However, the growth rate decreases as the thickness ratio increases and becomes more evident with the increase in middle layers. Furthermore, if the middle layer is composed of isotropic and anisotropic layers with hm and hc for each layer as shown in Figure 3(c), then the effect of different numbers of layers and thickness ratios (hm/hc) on the frequency parameters should be reconsidered, and the temperature rise is set from 0°C to 200°C. Figure 9(a) shows that when piezoelectric, anisotropic material, and isotropic material layers together form a symmetrical composite structure, the frequency characteristics differ from those of previous examples. For each order frequency, the thickness ratio has a positive effect on free vibration parameters, and the value is between 1 and 5. If the thickness ratio is less than 1 and greater than 5, the frequency will change in the opposite direction. In addition, the change trends in Figures 9(c) and 9(d) are similar with the rules in Case 1. The temperature rise still has a negative but negligible effect on the frequency parameters in this example.Figure 8
Variations of the frequency parameter Ω versus the thickness ratio (hm/hp) and the number of layers for an N-layer piezoelectric laminated beam (PZT-5A/Al/Al/…/PZT-5A).
(a)(b)Figure 9
Variations of the frequency parameter Ω versus the thickness ratio (hm/hc), number of layers, and temperature rise for an N-layer piezoelectric laminated beam (PZT-5A/material-D/Al/…/PZT-5A).
(a)(b)(c)(d)


## 5. Conclusion

This article mainly studied the free vibration of thermo-electro-elastic piezoelectric beams with the sandwich structure. To obtain more general characteristics and simplify the study of related parameters, we introduced the spring boundary to simulate the boundary force at both ends and integrated the phase and scattering relationships of the structure to achieve the total MRRM equation. By performing an algorithm based on the GSS for the MRRM equation, the natural frequency of the structure was calculated and compared with the FEM results and existing literature to verify the correctness of the proposed method. To further explore the effects of geometric factors and material parameters on the vibration of the structure, we carried out a series of analyzes of the geometric and material characteristics of the structure, especially the middle metal layer. These parameters included elastic modulus, the number of layers, thickness ratio, ply orientation of anisotropic material, different boundary conditions, each spring constant, temperature rise, and external electric potential of the beam. Through the coupling of external environmental parameters and geometric parameters, we arrived at the following conclusions:(1)
The free vibration parameters change obviously when the spring coefficient range is generally between 105 and 1010.(2)
For sandwich piezoelectric beams, the elastic modulus, thickness, and the number of layers of the middle metal layer have a positive effect on the free vibration.(3)
When the middle metal layer is composed of an anisotropic material, the ply orientation and lamination method influence the free vibration significantly. The more layers that change the ply orientation, the larger the range of frequency fluctuation.(4)
If the middle layer is composed of isotropic and anisotropic materials, then its thickness ratios have different effects on the frequency parameters in different numerical ranges, and additional layers will increase the frequency parameters.(5)
Factors in the external coupled environment, such as temperature and electric potential, have a negative effect on free vibration parameters of sandwich piezoelectric laminated beams. The temperature rise and external potential increase will cause a negligible decrease in frequency.



---

*Source: 1328886-2021-06-17.xml*