Coupled thermo-electro-elastic forced vibrations of piezoelectric laminated beams by means of Green's functions

International Journal of Mechanical Sciences 156 (2019) 355–369

Contents lists available at ScienceDirect

International Journal of Mechanical Sciences journal homepage: www.elsevier.com/locate/ijmecsci

Coupled thermo-electro-elastic forced vibrations of piezoelectric laminated beams by means of Green’s functions X. Zhao a,∗, F.J.N. Iegaink a, W.D. Zhu b,∗, Y.H. Li c a

School of Civil Engineering and architecture, Southwest Petroleum University, Chengdu 610500, PR China Department of Mechanical Engineering, University of Maryland, Baltimore County, Baltimore, MD 21250, USA c School of mechanics and engineering, Southwest Jiaotong University, Chengdu 610031, PR China b

a r t i c l e Keywords: Euler–Bernoulli beam Thermos-electro-elastic Energy harvesting Green’s function Coupling effects

i n f o

a b s t r a c t Increased interest in thermo-electro-elastic during the last decades can be assigned to the fact that the study of thermo-electro-elastic mechanical coupled behavior in smart structure. For instance, the coupled thermo-electroelastic dynamic factors are meaningful in designing process of some sensors and energy harvesting products. This paper presents the steady-state analytical solutions for the coupled thermo-electro-elastic forced vibrations of piezoelectric laminated beams. Two types of damping, material damping and air damping, are introduced to the coupled vibration system. Based on Erturk and Inman’ work, a novel electric circuit model is developed. Thus, for this model, the electric factor is introduced into the classical coupled thermoplastic vibration problems for the piezoelectric laminated beams. Using the decoupled method solve the three fields of coupling problems by developing a generalized form of Green’s functions. The interactions between the thermal and electric factors are obtained analytically. The convergence of the present solutions is firstly verified in the numerical section followed by the FEM results used to validate the achieved solutions. From the sample numerical calculations, it is seen that the two-dimensional temperature field presents different distributions for different electric conditions. Moreover, the optimum heat transfer coefficient is furnished for the energy harvesting problem of the laminated beam.

1. Introduction The artificial intelligent (AI) technology started out with dreams of human-level artificial general intelligence. During the last years, enthusiasm for this technology has risen in diverse engineering and scientific field with no exceptions. Concerning the structure material research area, several scholars have developed some intelligent materials and designed corresponding intelligent structures which can be used in many distributed sensors and actuators to make them have self-controlling and self-monitoring capabilities [1–4]. Thus, intelligent materials such as piezoelectric materials, electromagnetoelastic materials, shape memory alloys, electrostrictive materials, and pyro electric material, piezoelectric materials are applied more fashionably in many sensors and controllers of a wide range of mechanical systems and are utilized in medical sciences and aerospace structures, communications, printing industry and measurement technology. [5–8]. Before describing and discussing the vibration problems of piezoelectric structure, it’s worthwhile having a brief review of the main piezoelectric material property. One of the unique characteristics of this material is that it is reversible. It means that piezoelectric materials exhibit the direct piezoelectric effect which generates electricity when

∗

the external load is applied. On the contrary, the converse piezoelectric effect is the generation of stress when electric fields are applied. It should be noticed that in thermal environments, the electric potential difference can be also triggered by outside fluctuated temperature. This specific phenomenon is the so-called pyroelectric effect. [9–11]. This phenomenon changes thus the dynamic problems of piezoelectric structures from coupled electro-elastic vibrations to coupled thermo-electroelastic vibrations. At this stage, the thermal load must be considered and a tremendous amount of sensors are set in extremely complex environments. During the designing process of these piezoelectric structures, piezoceramics materials are mostly utilized. However, they are inorganic materials with low strength and bad quality factor. In order to enhance the strength and mechanical property, piezoceramics materials are merged with some rigid materials such as aluminum alloys (with higher strength and lightweight) to form some reliable structures. That is how the piezoelectric structures become composite structures or laminated structures constituted by several piezoceramics layers and subgrade layers. Erturk and Inman [12] studied the coupled vibration problems of the piezoelectric beams. By introducing a closed circuit for depicting the electric field in the piezoelectric beam, they developed the coupled

Corresponding author. E-mail addresses: [email protected] (X. Zhao), [email protected] (W.D. Zhu).

https://doi.org/10.1016/j.ijmecsci.2019.04.011 Received 11 November 2018; Received in revised form 14 March 2019; Accepted 5 April 2019 Available online 6 April 2019 0020-7403/© 2019 Elsevier Ltd. All rights reserved.

X. Zhao, F.J.N. Iegaink and W.D. Zhu et al.

International Journal of Mechanical Sciences 156 (2019) 355–369

electromechanical vibration equation of a piezoelectric beam. Subsequently, the model has been employed to investigate energy harvesting problems. Erturk and Inman [12] worked on the energy harvesting problems of the piezoelectric materials while Danesh-Yazdi et al. [13] investigated the energy harvesting problems of a unimorph Euler–Bernoulli piezoelectric harvester by means of Green’s function method, which is just suitable for a thin unimorph piezoelectric laminated beam. With a view to developing the limited electromechanical model, Dietl et al. [14] improved the Erturk and Inman’s model by using the Timoshenko assumption and studying the bimorph piezoelectric energy harvester, which is a thick beam structure with two piezoelectric layers. Recently, Zhao et al. [15] developed the Erturk and Inman’s model to the unimorph piezoelectric Timoshenko model by virtue of Green’s function method and obtained the closed-form solutions of the unimorph piezoelectric energy harvester. Guan and Liao [16] proposed a novel rotational piezoelectric energy harvester and studied it theoretically. Fan et al. [17] studied a type of multi-directional vibration piezoelectric energy harvester that obtain energy from both sway and bi-directional vibrations. Zhang et al. [18] investigated an efficiency rotational piezoelectric energy harvester for harvesting wind energy. With the assistant of experiments, Xiong and Wang [19] presented the applications of the piezoelectric energy harvesters in public roadways. The piezoelectric energy harvesters can scavenge the pavement deformation energy induced by the vehicles. To date investigations on the coupled vibration behavior in piezoelectric structures have mainly considered an interaction among electric and mechanical fields. However, investigations on thermal effect which should not be neglected in most situations have been very few. Mindlin [20] investigated the thermo-electro-mechanical dynamic problems of the piezoelectric plates and introduced the piezothermoelastic constitutive relation of the 3D linear thermopiezoelectric elastomer. Tzou and Bao [21] firstly proposed the thermo-electromechanical shell vibration equations by employing the piezothermoelastic constitutive relation of anisotropic piezoelectric materials and utilized the theory to study some sensing and controlling problems. Jiang and Li [22] proposed a novel finite element model to study the energy harvesting problems of the piezothermoelastic composite beams. Li et al. [23] investigated the effect of piezoelectric-thermo-elastic coupling on the piezoelectric material in piezoelectric energy harvester for the flexural vibration. In recent years, the research hotspots of the thermo-electromechanical dynamic problems and energy harvesting problems of the piezoelectric structures have switched to the functionally graded materials field and micro- or nanostructures field for obtaining more intelligent and more high-tech structures. Functionally graded materials (FGM) is specially designed for adapting the high temperature environments [24–26]. This typical characteristic makes FGM compatible with the piezothermoelastic structures which also work in the high temperature environments. Using the voltage of actuator to express the electric field, Komijani and Gracie [27] investigated the nonlinear thermo-electro-mechanical vibration behaviors of the Functionally Graded Piezoelectric Material (FGPM) beams under in-plane and out-of-plane thermal, mechanical and electric loads and studied the stable problems of the FGPM beams. Larkin and Abdelkefi [28] considered single-layered functionally graded piezoelectric energy harvesting systems and studied the performances of the FGM piezoelectric energy harvesters. Amini et al. [29] introduced finite element modeling methods to analyze FGM piezoelectric harvesters in unimorph or bimorph. Based on Euler-Bernoulli theory, Cao et al. [30] established the PFGM cantilever beam vibration models and optimized energy harvesting characteristics of PFGM cantilever beams. Resize of sensors and actuators are the consequence of the modern technologies development. Now they become smaller reaching micro or Nano level. Microscopically, the coupled multi-physical effects are amplified which lead the vibration of a micro piezoelectric beam to be commonly considered as multi-physical dynamic problems. Farzad and

Erfan [31] investigated the thermo-electro-mechanical dynamic behaviors of the FGM Nano beams under thermal load and applied voltage. Jung et al. [32] for his side designed a piezoelectric micro wind energy harvester to support the wireless sensor networks and optimized the energy output by experiments whereas Kuo et al. [33] created a bimorph piezoelectric cantilever micro energy generator and studied the energy harvesting for cell phone and some micro sensors. In brief, we have seen among the previous literature that most studies of the piezoelectric energy harvester just considered mechanical loads. In thermal environments, the external loads involved are thermal excitations. In this case, the vibrations of the piezoelectric energy harvesters should be piezoelectric-thermo-elastic coupled vibrations. To the best knowledge of authors, the piezoelectric beams subjected by thermal loads are very limited in the energy harvesting research. No one has yet attempt to find the optimized heat transfer coefficient 𝛼 1 to make the voltage larger. Moreover, the analytical solutions of the coupled thermoelectro-elastic forced vibrations of piezoelectric laminated beams with two damping effects have not yet been reported in the literature. In this paper, the coupled thermo-electro-elastic dynamic problems and the energy harvesting studies of a laminated piezoelectric beam with a piezoceramic layer and a metallic substructure layer is investigated. The intention is to develop the steady-state responses of the coupled thermo-electro-elastic forced vibrations of piezoelectric beams subjected to an external thermal load. Only the temperature factor in the metallic layer (Aluminum alloys) will be considered due to the fact that its thermal conductivity is bigger that piezoceramics materiel. The temperature distribution of the sublayer is presented by a two-dimension heat transfer equation, and it is the classical Euler–Bernoulli beam theory that is applied to build the coupled vibration equation. A developed electromechanical model based on the Erturk and Inman’s work [12] is developed to illustrate the electric field of the beam. According to the work of Zhao et al. [34, 35], the explicit expressions of the displacement, temperature and voltage can be presented. In the numerical section, the present solutions are verified by comparing the archived result with the FEM results. The effects of the convective heat transfer coefficient to the energy harvesting and the influences of the electric resistance on the temperatures are chiefly discussed. Additionally, the effects of heating positions, material damping and air damping are also investigated. This study provides some original thoughts to help design some novel electrical sensors and temperature sensors and gives some optimizing strategies to small energy harvesting products.

2. Governing equations 2.1. Coupled vibration equation for the piezoelectric laminated beam Fig. 1(a) exhibits a simple supported unimorph piezoelectric laminated beam subjected to a harmonic heat load q (x, y, t). The piezoelectric laminated beam, which is considered to be a uniform Euler– Bernoulli beam, is fabricated by the upper piezoceramic PZT 5A & 5H layer and the lower substructure metal layer. The metal layer is generally made by Aluminum alloys which have an advantage of high strength and light weight. It is assumed that the piezoelectric and substructure layers are perfectly bonded. The length of the cantilever beam is L, and the R l denotes the electric resistance. The electric circuit shown in Fig. 1(a) is attached to the piezoelectric layer. The cross-section of the piezoelectric laminated beam is shown in Fig. 1(b). In the Fig. 1(b), the symbol bs (bp ) denotes the width of the metal sublayer (PZT 5A & 5H layer); the hp (hs ) is the thickness of the PZT 5A & 5H layer (metal sublayer); the total height of the beam is h = hp + hs ; the size ha stands for the distance from the neutral axis (N. A.) to the lowest surface of the metallic sublayer, and the hb is the distance from N. A. to the bottom of the PZT 5A & 5H layer; the distance from N. A. to the top surface of the PZT 5A & 5H layer is expressed by symbol hc . 356

X. Zhao, F.J.N. Iegaink and W.D. Zhu et al.

International Journal of Mechanical Sciences 156 (2019) 355–369

Fig. 1. A simple supported piezoelectric laminated beam subjected to a point heat load (a) and cross section (b).

initial temperature. Moreover, the effect of the temperature distribution of the interface between PZT and metallic layer 𝑇̄ (𝑥, ℎ𝑏 , 𝑡) is considered in the stress-strain relation of the PZT layer. Based on Eq. (2), the expression of the bending moment M (x, t) can be written as

The heat transfer capacity difference between the two materials (1.5 W m− 1 K−1 [36]) and (200–237 W m− 1 K−1 ) respectively for the PZT 5A & 5H and the aluminum alloys along with and their thickness (thinner in PZT layer) lead to the following reasonable assumptions which help us simplify the coupled dynamic problems: (1) The heat transfer process of the PZT 5A &5H layer is not considered in this study; (2) The interface of the PZT layer and the metal layer is seen as heat insulation; These assumptions mean that thermal conduction occurs only in the metal layer and we only consider the electrical effect in the PZT layer. The metal layer loaded by a harmonic thermal flux will produce a thermal stress. Then, the thermal deformation of the metal layer induces the deflection of the PZT layer. Based on the piezoelectric effects, the voltage will occur in the circuit Based on the Euler Bernoulli beam model and the Ref. [12], the governing equation of the laminated beam can be written as 𝜕2 𝑀 𝜕 𝑥2

+ 𝑐𝑠 𝐼

𝜕5 𝑤 𝜕 𝑥4 𝜕𝑡

+ 𝑐𝑎

𝜕2 𝑤

𝜕𝑤 +𝑚 =0 𝜕𝑡 𝜕 𝑡2

𝑀(𝑥, 𝑡) = −

𝑠 𝜎𝑥𝑥 𝑏𝑠 𝑦𝑑𝑦 −

ℎ𝑐

∫ℎ𝑏

𝑝 𝜎𝑥𝑥 𝑏𝑝 𝑦𝑑𝑦.

(3)

Substituting the constitutive relation Eq. (2) into Eq. (3) obtains 𝑀(𝑥, 𝑡) =

ℎ𝑏

∫ℎ𝑎

𝑏𝑠 𝐸𝑠 𝑦2 ℎ𝑏

+

∫ℎ𝑎

ℎ

ℎ

𝑐 𝑐 𝑣(𝑡) 𝜕2 𝑤 𝜕2 𝑤 𝑑𝑦 + 𝐸𝑝 𝑏𝑝 𝑦2 − 𝑑31 𝐸𝑝 𝑏𝑝 𝑦 𝑑𝑦 ∫ℎ𝑏 ℎ𝑝 𝜕 𝑥2 𝜕 𝑥2 ∫ℎ𝑏

𝛼𝑇 𝐸𝑠 𝑏𝑠 𝑦(𝑇̄ − 𝑇0 )𝑑𝑦 +

ℎ𝑐

∫ℎ𝑏

𝛼𝑇 𝐸𝑝 𝑏𝑝 𝑦𝑇̄ (𝑥, ℎ𝑏 , 𝑡)𝑑𝑦

(4)

where the electric field strength Ey =−v (t)/hp, and the v (t) is the voltage of the PZT layer. The equation Eq. (4) is still quite complex and after some manipulations, a more useful form can be rewritten as

(1)

where the w (x, t) is the transversal displacement; the M (x, t) expresses the bending moment; the symbol m denotes the mass per unite length of the beam; the c s and c a respectively stand for the inner material damping and outer air damping of the dynamic system. For the aim of obtaining the expression of the bending moment M (x, t), the following constitutive relation is taken into consideration 𝑝 𝑠 𝜎𝑥𝑥 = 𝐸𝑝 [𝜀𝑝𝑥𝑥 − 𝑑31 𝐸𝑦 − 𝛼𝑇 𝑇̄ (𝑥, ℎ𝑏 , 𝑡)], 𝜎𝑥𝑥 = 𝐸𝑠 [𝜀𝑠𝑥𝑥 − 𝛼𝑇 (𝑇̄ − 𝑇0 )],

ℎ𝑏

∫ℎ𝑎

𝑀(𝑥, 𝑡) = (𝐸𝐼)𝑒𝑓 𝑓

𝜕2 𝑤 + 𝜗𝑣(𝑡)[𝐻(𝑥 − 𝑥1 ) − 𝐻(𝑥 − 𝑥2 )] 𝜕 𝑥2 ℎ

+

𝑏 1 2 (ℎ − ℎ2𝑏 )𝛼𝑇 𝐸𝑝 𝑏𝑝 𝑇̄ (𝑥, ℎ𝑏 , 𝑡) + 𝛼 𝐸 𝑏 𝑦(𝑇̄ − 𝑇0 )𝑑𝑦 ∫ℎ𝑎 𝑇 𝑠 𝑠 2 𝑐

(5)

where the expressions of the effective transverse stiffness (EI) eff and the coupling coefficient ϑ are represented as

(2)

where the 𝜎 p xx (𝜀 p xx) stands for normal stress (normal strain) of the PZT layer and 𝜎 s xx (𝜀 s xx) is normal stress (normal strain) of the metallic layer; the Ep and Es are Young’s modulus of the PZT layer and the metal layer, respectively; the d31 denotes piezoelectric constants; and the y-axial electric field is represented by the symbol Ey ; the 𝛼 T illustrates heat transfer coefficient of the metallic layer; the 𝑇̄ (𝑥, 𝑦, 𝑡) expresses the temperature distribution of the metal layer, and the T0 is the

(𝐸𝐼)𝑒𝑓 𝑓 =

𝑏𝑠 𝐸𝑠 (ℎ3𝑏 − ℎ3𝑎 ) + 𝑏𝑝 𝐸𝑝 (ℎ3𝑐 − ℎ3𝑏 ) 3

,𝜗=−

𝐸𝑝 𝑑31 𝑏𝑝 2ℎ𝑝

(ℎ2𝑐 − ℎ2𝑏 ).

(6)

The interval range x1 ≤ x ≤ x2 denotes the covered region of the piezoelectric layer. For example, the interval 0 ≤ x ≤ L means that the piezoelectric layer covers the entire length of the beam. 357

X. Zhao, F.J.N. Iegaink and W.D. Zhu et al.

International Journal of Mechanical Sciences 156 (2019) 355–369

Inserting Eq. (5) into Eq. (1), we can derive the forced vibration equation of the piezoelectric laminated beam as 𝜕4 𝑤 𝜕5 𝑤 𝜕𝑤 𝜕2 𝑤 + 𝑐𝑠 𝐼 + 𝑐𝑎 +𝑚 𝜕𝑡 𝜕 𝑥4 𝜕 𝑥4 𝜕𝑡 𝜕 𝑡2 [ ] ℎ𝑏 d𝛿(𝑥 − 𝑥1 ) d𝛿(𝑥 − 𝑥2 ) 𝜕 2 𝑇̄ +𝜗 ⋅ 𝑣(𝑡) − =− 𝛼 𝐸𝑏 𝑦d𝑦 ∫ℎ𝑎 𝑇 𝑠 𝑠 𝜕 𝑥2 d𝑥 d𝑥 1 − (ℎ2𝑐 − ℎ2𝑏 )𝛼𝑇 𝐸𝑝 𝑏𝑝 𝑇̄ (𝑥, ℎ𝑏 , 𝑡). 2

𝜕4 𝑤 𝜕5 𝑤 𝜕𝑤 + 𝑐𝑠 𝐼 + 𝑐𝑎 𝜕𝑡 𝜕 𝑥4 𝜕 𝑥4 𝜕𝑡 [ ] ℎ𝑏 d𝛿(𝑥 − 𝑥1 ) d𝛿(𝑥 − 𝑥2 ) 𝜕2 𝑤 𝜕 2 𝑇̄ +𝑚 + 𝜗 ⋅ 𝑣(𝑡) − =− 𝛼𝑇 𝐸𝑠 𝑏𝑠 𝑦d𝑦 2 ∫ℎ𝑎 d𝑥 d𝑥 𝜕𝑡 𝜕 𝑥2 1 − (ℎ2𝑐 − ℎ2𝑏 )𝛼𝑇 𝐸𝑝 𝑏𝑝 𝑇̄ (𝑥, ℎ𝑏 , 𝑡), 2 (13a)

(𝐸𝐼)𝑒𝑓 𝑓

(𝐸𝐼)𝑒𝑓 𝑓

(7)

( 2.2. Equation for electric circuit in PZT layer

𝜆𝑇

Since the temperature distribution of the interface between the PZT and metal layer influences the constitution model of the PZT layer, the constitutive relation of the PZT layer is different from that of Erturk and Inman [12]. The model can be written as 𝐷𝑦 = 𝑑31 𝐸𝑝 𝜀𝑥𝑥 − 𝜀𝑠33

𝑣(𝑡) + 𝑝𝑇 𝑇̄ (𝑥, ℎ𝑏 , 𝑡) ℎ𝑝

𝐶𝑝

𝜀𝑠33 𝑏𝑝 (𝑥2 − 𝑥1 ) ℎ𝑝

, 𝛽 = 𝐸𝑝 𝑑31 𝑏𝑝 ℎ𝑝𝑐 = −

2ℎ𝑝 ℎ𝑝𝑐 ℎ2𝑐 − ℎ2𝑏

𝜗.

(14)

𝑊 ′ + 𝐶𝐷1 𝑊 + 𝐶𝐷2 ⋅ 𝑉 [𝛿 ′ (𝑥 − 𝑥1 ) − 𝛿 ′ (𝑥 − 𝑥2 )] = 𝐶𝑇 1

ℎ𝑏

∫ℎ𝑎

𝑇 ′′ 𝑦d𝑦 + 𝐶𝑇 2 𝑇 (𝑥, ℎ𝑏 ),

Δ𝑇 + iΩ𝐶𝑇 3 𝑇 = −𝐶𝐷4 𝑊 ′′ − 𝜆−1 𝑇 𝑞 (𝑥, 𝑦),

(10)

𝑉 = 𝐶𝐷3

where 𝐶𝑝 =

(13c)

where the W (x), V and T (x, y) are respectively the steady-state transversal deflection, steady-state voltage and steady-state two-dimension temperature. Substituting Eq. (14) into Eq. (13), we arrive at

(9)

𝑥2 3 𝑥2 𝑑𝑣(𝑡) 𝑣(𝑡) 𝜕 𝑤(𝑥, 𝑡) + = −𝛽 𝑑𝑥 + 𝑝 𝑏 𝑇̄ (𝑥, ℎ𝑏 , 𝑡)𝑑𝑥. 2 ∫𝑥1 𝑇 𝑝 ∫𝑥1 𝑑𝑡 𝑅𝑙 𝜕 𝑥 𝜕𝑡

𝑥2 3 𝑥2 𝑑𝑣(𝑡) 𝑣(𝑡) 𝜕 𝑤(𝑥, 𝑡) + = −𝛽 𝑑𝑥 + 𝑝 𝑏 𝑇̄ (𝑥, ℎ𝑏 , 𝑡)𝑑𝑥. ∫𝑥1 𝑇 𝑝 ∫𝑥1 𝑑𝑡 𝑅𝑙 𝜕 𝑥2 𝜕𝑡

𝑤(𝑥, 𝑡) = 𝑊 (𝑥)𝑒i 𝜔𝑡 , 𝑣(𝑡) = 𝑉 𝑒i 𝜔𝑡 , 𝑇̄ (𝑥, 𝑦, 𝑡) = 𝑇 (𝑥, 𝑦)𝑒i 𝜔𝑡

In the Eq. (9), the h pc denotes the distance from N. A. to the central plane of the piezoelectric layer, which is shown in Fig. 1(b). Utilizing the constitutive relation Eq. (8) and the method in Ref. [12], the electric circuit equation with displacement coupling can be given as 𝐶𝑝

) ( ) 𝛼𝑇2 𝐸𝑠 𝑇0 𝜕 𝑇̄ 𝛼𝑇 𝐸𝑠 𝑇0 ℎ𝑝𝑐 𝜕 𝜕 2 𝑤 𝜕 2 𝑇̄ 𝜕 2 𝑇̄ 1 = + + + 𝑞 (𝑥, 𝑦, 𝑡), 𝜆𝑇 𝜕𝑡 𝜆𝑇 𝜕𝑡 𝜕 𝑥2 𝜆𝑇 𝜕 𝑥2 𝜕 𝑦2

Assuming the heat load is harmonic i.e. q (x, y, t) = q (x, y) ei𝜔 t , the corresponding transversal displacement, voltage and temperature distribution can be given the same assumption

(8)

𝜕 2 𝑤(𝑥, 𝑡) . 𝜕 𝑥2

+

(13b)

where the Dy (x, t) is electric displacement along the thickness of the beam, and the 𝜀 s 33 is dielectric constant. Moreover, the pT is the pyroelectric constants. The function 𝜀xx (x, t) stands for the average bending strain. On the basis of Euler–Bernoulli assumption, the 𝜀xx (x, t) should be assumed as [12] 𝜀𝑥𝑥 = −ℎ𝑝𝑐

𝑐𝑝

𝑥2

∫𝑥1

𝑊 ′′ 𝑑𝑥 +

𝑥2

∫𝑥1

𝑝𝑇 𝑏𝑝 𝑇 (𝑥, ℎ𝑏 )𝑑𝑥.

(15a) (15b) (15c)

In Eq. (15),

(11)

𝐶𝐷1 =

2.3. Heat conduction equation for metallic layer

i𝜔𝑐 𝑎 − 𝜔2 𝑚 i𝜔𝜗 ,𝐶 = , (𝐸𝐼)𝑒𝑓 𝑓 + i𝜔𝑐𝑠 𝐼 𝐷2 (𝐸𝐼)𝑒𝑓 𝑓 + i𝜔𝑐𝑠 𝐼

−𝑐𝑝 − (𝛼𝑇 )2 𝐸𝑠 𝑇0 −𝐸𝑠 𝑏𝑠 𝛼𝑇 1 , 𝐶 = − (ℎ2𝑐 − ℎ2𝑏 )𝛼𝑇 𝐸𝑝 𝑏𝑝 , 𝐶𝑇 3 = , (𝐸𝐼)𝑒𝑓 𝑓 + i𝜔𝑐𝑠 𝐼 𝑇 2 2 𝜆𝑇 i𝜔𝛼𝑇 𝐸𝑠 𝑇0 ℎ𝑝𝑐 −i𝜔𝛽𝑅𝑙 = ,𝐶 = . i𝜔𝐶𝑝 𝑅𝑙 + 1 𝐷4 𝜆𝑇

𝐶𝑇 1 =

It should be clear that the metallic layer of the piezoelectric laminated beam has no ability to produce electric, there is only a thermal factor coupled with the structure vibration in this layer. Further insights into the coupled thermoelastic vibrations of some beams have been studied in some references, and the governing equations of the coupled thermoelastic vibrations of the beams have been established in these papers [37–39]. According to these Refs. [37–39], the heat transfer equation with displacement coupling effect can be obtained as (

𝐶𝐷3

(16) Additionally, in Eq. (15), the prime indicates a derivative with respect to the spatial coordinate x, and the symbol ∆, which is equal to 𝜕 2 /𝜕 x2 + 𝜕 2 /𝜕 y2 , denotes the Laplace operator. The sign “i” stands for the imaginary unit. If the piezoelectric effect is ignored i.e. the upper layer is not PZT material but the same metallic material with the sublayer, d31 = 𝜀 s 33 = 0, the Eq. (15) can be degenerated to the coupled thermoelastic vibration equations of the Euler–Bernoulli beams [35,40]. Furthermore, if the thermal effect is not considered, the Eq. (15) will be reduced to the coupled electromechanical vibration equations of the Euler–Bernoulli beams [12].

) ( ) 𝛼𝑇2 𝐸𝑠 𝑇0 𝜕 𝑇̄ 𝛼𝑇 𝐸𝑠 𝑇0 ℎ𝑝𝑐 𝜕 𝜕 2 𝑤 𝜕 2 𝑇̄ 𝜕 2 𝑇̄ 1 + = + + + 𝑞 (𝑥, 𝑦, 𝑡) 𝜆𝑇 𝜆𝑇 𝜕𝑡 𝜆𝑇 𝜕𝑡 𝜕 𝑥2 𝜆𝑇 𝜕 𝑥2 𝜕 𝑦2 𝑐𝑝

(12) where the symbol 𝜆T , 𝛼 T , cp are respectively the coefficients of thermal conductivity, the thermal expansion and the heat capacity per unit volume. The function q (x, y, t) is the external heat load, and the range of y is h a ≤ y ≤ hb . Moreover, in the Eq. (12), the term (𝛼 T Es T0 hpc /𝜆T ) · 𝜕 (𝜕 2 w/𝜕 x2 )/𝜕 tis the displacement coupling term.

3. Green’s functions for the steady-state heat transfer equations The Green’s functions of the heat transfer equations with various boundary conditions can be obtained according Ref [35,40] with this approach of eigenfunction expansion based on Sturm–Liouville theory [41]. For seek of saving space, the Green’s functions will just be introduced. Nevertheless, the detailed solving process can be seen in the work of Zhao et al. [35].

2.4. Coupled thermo-electro-mechanical governing equations Reorganizing the Eqs. (7), (10) and (12) obtains the coupled thermos-electro-mechanical vibration equations of the piezoelectric laminated beam 358

X. Zhao, F.J.N. Iegaink and W.D. Zhu et al.

International Journal of Mechanical Sciences 156 (2019) 355–369

The expression of the Green’s function of the heat conduction equation Eq. (15b) can be written as the following form [35] 𝐺𝑇 (𝑥, 𝑦; 𝑥0 , 𝑦0 ) =

Utilizing the Laplace transform method, we can derive the Green’s functions

( ) ( ) ( ) ∞ ∑ ∞ ∑ 2 sin 𝑛𝜋𝑥0 ∕𝐿 Φ𝑚 𝑘𝑚 , 𝑦0 sin (𝑛𝜋𝑥∕𝐿)Φ𝑚 𝑘𝑚 , 𝑦 ( ) 𝑞𝑚 𝐿 𝑘2𝑚 + 𝑛2 𝜋 2 ∕𝐿2 − i𝜔𝐶𝑇 3 𝑚=1 𝑛=1

𝐺𝐷1 (𝑥; 𝑥0 ) = 𝐻(𝑥 − 𝑥0 )𝜙11 (𝑥 − 𝑥0 ) + 𝜙2 (𝑥)𝑊 (0) + 𝜙3 (𝑥)𝑊 ′ (0) + 𝜙4 (𝑥)𝑊 ′′ (0) + 𝜙5 (𝑥)𝑊 ′′′ (0),

(23)

(17) and

where the eigenfunction Φ (∙, ∙), eigenvalue k m and q m are determined by the boundary conditions of the two-dimension heat conductive region. For different boundary conditions, the expressions of the Φ (∙, ∙), k m and q m have been listed in the Ref. [35]. Now, a heat load q (x, y, t) = q0 𝛿 (x – x0 ) 𝛿 (y – y0 ) ei𝜔 t is taken into account. Based on the Green’s function Eq. (17) and superposition principle, the temperature solution of the heat transfer equation Eq. (15b) can be illustrated as

𝐺𝐷2 (𝑥; 𝑥0 ) = 𝐻(𝑥 − 𝑥0 )𝜙12 (𝑥 − 𝑥0 ) + 𝜙2 (𝑥)𝑊 (0) + 𝜙3 (𝑥)𝑊 ′ (0) + 𝜙4 (𝑥)𝑊 ′′ (0) + 𝜙5 (𝑥)𝑊 ′′′ (0),

where the H (∙) is unit step function. The unknown constants W (0), W′ (0), W″ (0) and W‴ (0) can be determined by the boundary conditions of the beam. To be more specific terms,

( ) ( ) 𝑇 𝑥, 𝑦; 𝑥0 , 𝑦0 = 𝑇1 (𝑥, 𝑦) + 𝜆−1 𝑇 𝑞0 𝑇2 𝑥, 𝑦; 𝑥0 , 𝑦0 , 𝑇1 (𝑥, 𝑦) =

∞ ∑ ∞ ∑ 𝑚=1 𝑛=1

( ) 𝐹01 (𝑚, 𝑛) sin (𝑛𝜋𝑥∕𝐿)Φ𝑚 𝑘𝑚 , 𝑦

𝐿

∫0

(24)

𝜙11 (𝑥) =

4 ∑ 𝑖=1

𝑊 ′′ (𝜉) sin (𝑛𝜋𝜉∕𝐿)d𝜉, 𝜙3 (𝑥) =

∞ ∑ ∞ ( ) ∑ ( ) ( ) 𝑇2 𝑥, 𝑦; 𝑥0 , 𝑦0 = 𝐹02 𝑚, 𝑛; 𝑥0 , 𝑦0 sin (𝑛𝜋𝑥∕𝐿)Φ𝑚 𝑘𝑚 , 𝑦

4 ∑ 𝑖=1

𝐴𝑖 (𝑥), 𝜙12 (𝑥) =

𝐴𝑖 (𝑥)𝑠2𝑖 , 𝜙4 (𝑥) =

4 ∑

𝐴𝑖 (𝑥) ⋅ (−𝑠𝑖 ), 𝜙2 (𝑥) =

𝑖=1 4 ∑ 𝑖=1

𝐴𝑖 (𝑥)𝑠𝑖 , 𝜙5 (𝑥) =

4 ∑ 𝑖=1

4 ∑ 𝑖=1

𝐴𝑖 (𝑥)𝑠3𝑖 , (25)

𝐴𝑖 (𝑥).

𝑚=1 𝑛=1

In the Eq. (25),

(18) where

) ( ℎ 2𝐶𝐷4 ∫ℎ 𝑏 Φ𝑚 𝑘𝑚 , 𝜂 d𝜂 𝑎 1 𝐹0 (𝑚, 𝑛) = ( ), 𝑞𝑚 𝐿 𝑘2𝑚 + 𝑛2 𝜋 2 ∕𝐿2 − i𝜔𝐶𝑇 3 ( ) ( ) 2 sin 𝑛𝜋𝑥0 ∕𝐿 Φ𝑚 𝑘𝑚 , 𝑦0 ( ) 𝐹02 𝑚, 𝑛; 𝑥0 , 𝑦0 = ( ). 𝑞𝑚 𝐿 𝑘2𝑚 + 𝑛2 𝜋 2 ∕𝐿2 − i𝜔𝐶𝑇 3

𝐴1 (𝑥) = 𝐴3 (𝑥) =

(19)

𝑊 (𝑥) = 𝑉 ⋅ 𝐶𝐷2 ⋅ 𝐺𝐷𝑖𝑠 (𝑥) + 𝐶𝑇 1

3.1. Green’s functions for the steady-state forced vibration equations with two damping effects

+ 𝐶𝑇 2

𝑊 + 𝐶𝐷1 𝑊 = 𝐶𝑇 1

∫ℎ𝑎

𝑇 𝑦d𝑦 + 𝐶𝑇 2 𝑇 (𝑥, ℎ𝑏 ) − 𝐶𝐷2

(21)

𝐺𝐷1 (𝑥; 𝜉)𝑇 ′′ (𝜉, 𝑦)𝑦d𝑦𝑑𝜉

𝐺𝐷1 (𝑥; 𝜉)𝑇 (𝜉, ℎ𝑏 )𝑑𝜉.

(27)

[

and 𝑊 ′ 𝑐 + 𝐶𝐷1 𝑊 = 𝛿 ′ (𝑥 − 𝑥0 ).

∫0 ∫ℎ𝑎

In this section, we are about to decoupling of the coupled thermoselectro-elastic dynamic problem. To undertaking this, the coupled dynamic system Eq. (15) will be decoupled to derive the explicit expressions of the steady-state displacement, the voltage, and the temperature distribution. Remember that we have seen in Section 3 the heat load of the system is the point heat load q (x, y) = q0 𝛿 (x – x0 ) 𝛿 (y – y0 ) e i𝜔 t , and the length of the PZT layer is assumed to be equal to L. Initially, according to the Eqs. (15c), (18) and (27), the expression of the voltage can be provided

(20)

It can be shown in the Eq. (20) that there are two kind of load ℎ which are the heat load 𝐶𝑇 1 ∫ℎ 𝑏 𝑇 ′′ 𝑦d𝑦and the electric load CD 2 V [𝛿′ 𝑎 (x – x1 )−𝛿′ (x – x2 )]. Therefore, for the intent of obtaining the solutions of the Eq. (20), two Green’s functions GD 1 (x; x0 ) and GD 2 (x; x0 ) are defined, and according to the physical meaning of the Green’s functions [42], the GD 1 (x; x0 ) and GD 2 (x; x0 ) are respectively the solutions of the following equations 𝑊 ′ 𝑐 + 𝐶𝐷1 𝑊 = 𝛿(𝑥 − 𝑥0 ),

ℎ𝑏

4. Decoupling of the coupled thermos-electro-elastic dynamic problem

′′

⋅ 𝑉 [𝛿 ′ (𝑥 − 𝑥1 ) − 𝛿 ′ (𝑥 − 𝑥2 )].

𝐿

∫0

𝐿

where GDis (x) = GD 2 (x; x2 )−GD 2 (x; x1 ).

The coupled steady-state forced vibration equation of the Euler– Bernoulli beam with two damping has been derived and shown in the Eq. (15a). Reorganizing the Eq. (15a) can be obtained ℎ𝑏

e𝑠3 𝑥 e𝑠4 𝑥 , 𝐴 (𝑥) = . (𝑠3 − 𝑠1 )(𝑠3 − 𝑠2 )(𝑠3 − 𝑠4 ) 4 (𝑠4 − 𝑠1 )(𝑠4 − 𝑠2 )(𝑠4 − 𝑠3 ) (26)

In the Eqs. (25) and (26), the s i (i = 1, 2, 3, 4) are four characteristic roots of the characteristic equation s4 + CD1 = 0. The specific determination of the Green’s functions of various boundary conditions can be found in Ref. [43]. On the basis of the superposition principle, the steady-state displacement expression of the laminated beam can be obtained

As mentioned above, a more detailed process can be seen in Ref. [35]. In the Eq. (18), the temperature field T (x, y; x0 , y0 ) is constituted by two parts T1 (x, y) and 𝜆−1 Tq0 T2 (x, y; x0 , y0 ). The temperature field T1 (x, y) is produced by the coupled displacement heat load CD 4 y W″ (x). Physically, it means that the T1 (x, y) stands for the displacement coupling effect. Moreover, the temperature distribution T2 (x, y; x0 , y0 ) is the same as the temperature field of the uncoupled heat conduction equation due to the fact it is induced by external heat load.

′

e𝑠1 𝑥 e𝑠2 𝑥 , 𝐴 (𝑥) = , (𝑠1 − 𝑠2 )(𝑠1 − 𝑠3 )(𝑠1 − 𝑠4 ) 2 (𝑠2 − 𝑠1 )(𝑠2 − 𝑠3 )(𝑠2 − 𝑠4 )

(22)

𝑉 =

Referring to the work of Li et al. [43], the solutions of the Eqs. (21) and (22) can be obtained by means of Laplace transform. Because the specific deducing process have been presented in Ref. [43], for saving space, it is no need to list the process again.

𝐿∑ 𝐿∑ ′′ ′′ 𝐶𝑇 1 𝐶𝐷3 ∫0 2 (𝑚, 𝑛; 𝜉)𝑊 (𝜉)𝑑𝜉 + 𝐶𝑇 2 𝐶𝐷3 ∫0 4 (𝑚, 𝑛; 𝜉)𝑊 (𝜉)𝑑𝜉 ] 𝐿∑ ′′ + 𝑝 𝑇 𝑏 𝑝 ∫0 6 (𝑚, 𝑛; 𝜉)𝑊 (𝜉)𝑑𝜉

+

𝜆−1 𝑞 𝑇 0

[

𝐶𝑇 1 𝐶𝐷3

∑

1 − 𝐶𝐷3 ⋅ 𝐶𝐷2 ⋅ 𝐶𝐺1 ] ∑ ∑ 1 (𝑚, 𝑛; 𝑥0 , 𝑦0 ) + 𝐶𝑇 2 𝐶𝐷3 3 (𝑚, 𝑛; 𝑥0 , 𝑦0 ) + 𝑝𝑇 𝑏𝑝 5 (𝑚, 𝑛; 𝑥0 , 𝑦0 ) 1 − 𝐶𝐷3 ⋅ 𝐶𝐷2 ⋅ 𝐶𝐺1

(28) 359

X. Zhao, F.J.N. Iegaink and W.D. Zhu et al.

International Journal of Mechanical Sciences 156 (2019) 355–369 8

Inserting Eqs. (18) and (28) into Eq. (27) give rise to

Present solution FEM solution

𝑊 (𝑥) = ∫0

+ 𝐶𝑜𝑒𝑓3 (𝑥) + 𝐶𝑇 2

𝐿

∫0

∑ (𝑚, 𝑛; 𝜉)𝑊 ′′ (𝜉)d𝜉 + 𝐶 𝑜𝑒𝑓2 (𝑥) 2 𝐿∑

∫0

(𝑚, 𝑛; 𝜉)𝑊 ′′ (𝜉)d𝜉+𝐶𝑇 1

6

𝐿

∫0

𝐿

∫0

∑ (𝑚, 𝑛; 𝜉)𝑊 ′′ (𝜉)d𝜉

6

4

W(X;0.5,0.5)×109

𝐶 𝑜𝑒𝑓1 (𝑥)

𝐿

∑ (𝑚, 𝑛; 𝑥, 𝜉)𝑊 ′′ (𝜉)d𝜉 8

∑ (𝑚, 𝑛; 𝑥, 𝜉)𝑊 ′′ (𝜉)d𝜉 + 𝑁𝑜𝑛(𝑥; 𝑥0 , 𝑦0 )

(29)

10

2

which is actually a Fredholm integral equation [44]. Thus, the integral equation Eq. (29) actually realizes the decoupling of the original coupled differential system Eq. (15). In the Eqs. (28) and (29), [ ∑ 𝜆−1 𝑞 𝐶 ⋅ 𝐺 (𝑥) 𝐶𝑇 1 𝐶𝐷3 1 (𝑚, 𝑛; 𝑥0 , 𝑦0 ) 𝑇 0 𝐷2 ∑ 𝐷𝑖𝑠 ] ∑ +𝐶𝑇 2 𝐶𝐷3 3 (𝑚, 𝑛; 𝑥0 , 𝑦0 ) + 𝑝𝑇 𝑏𝑝 5 (𝑚, 𝑛; 𝑥0 , 𝑦0 ) 𝑁𝑜𝑛(𝑥; 𝑥0 , 𝑦0 ) = , 1 − 𝐶𝐷3 ⋅ 𝐶𝐷2 ⋅ 𝐶𝐺1 ∑ ∑ + 𝐶𝑇 1 𝜆−1 (𝑚, 𝑛; 𝑥0 , 𝑦0 , 𝑥) + 𝐶𝑇 2 𝜆−1 (𝑚, 𝑛; 𝑥0 , 𝑦0 , 𝑥) 𝑇 𝑞0 𝑇 𝑞0 7

4

0 0.00

0.25

0.50

0.75

1.00

X

(a)

9

0.50

𝐶 𝐶 𝐶 𝐺 (𝑥) 𝐶 𝐶 𝐶 ⋅ 𝐺𝐷𝑖𝑠 (𝑥) 𝐶𝑜𝑒𝑓1 (𝑥) = 𝐷3 𝑇 1 𝐷2 𝐷𝑖𝑠 , 𝐶𝑜𝑒𝑓2 (𝑥) = 𝑇 2 𝐷3 𝐷2 , 1 − 𝐶𝐷3 𝐶𝐷2 𝐶𝐺1 1 − 𝐶𝐷3 ⋅ 𝐶𝐷2 ⋅ 𝐶𝐺1 𝑝𝑇 𝑏𝑝 𝐶𝐷2 ⋅ 𝐺𝐷𝑖𝑠 (𝑥) 𝐶𝑜𝑒𝑓3 (𝑥) = . (30) 1 − 𝐶𝐷3 ⋅ 𝐶𝐷2 ⋅ 𝐶𝐺1

-0.002500 0.01819 0.03888 0.05956 0.08025 0.1009 0.1216 0.1423 0.1630

0.25

𝐿

Y

where 𝐶𝐺1 = ∫0 [𝐺𝐷2 ′′ (𝑥; 𝐿) − 𝐺𝐷2 ′′ (𝑥; 0)]𝑑𝑥. The expressions of the Σi (i = 1, 2 … 10) can be seen in Appendix A. It is assumed that the orders of the double series in Eq. (29) are the finite numbers r1 and r2 which are large enough to ensure the solutions of Eq. (29) convergent. Then the integral equation Eq. (29) can be recasted to the following form ) ∞ ∑ ∞ ( ∑ 𝑛2 𝜋 2 𝑊 (𝑥; 𝑥0 , 𝑦0 ) = 𝐶𝑇 1 − 𝐹01 (𝑚, 𝑛)𝑓𝐷1 (𝑛, 𝑥)𝐹11 (𝑚) 𝐿2 𝑚=1 𝑛=1

-0.25

-0.50 0.00

𝐿

× ∫0 𝑊 ′′ (𝜉) sin (𝑛𝜋𝜉∕𝐿)d𝜉 ) ∞ ∑ ∞ ( ∑ 𝑛2 𝜋 2 + 𝐶𝑜𝑒𝑓1 (𝑥) − 𝐹01 (𝑚, 𝑛)𝑓𝐷 (𝑛)𝐹11 (𝑚) 𝐿2 𝑚=1 𝑛=1 ×

𝐿 ∫0

𝑊

+ 𝐶𝑜𝑒𝑓2 (𝑥) 𝐿

× ∫0 𝑊

𝐹01 (𝑚, 𝑛)𝑓𝐷 (𝑛)Φ𝑚 𝑚=1 𝑛=1 ′′ (𝜉) sin (𝑛𝜋𝜉∕𝐿)d𝜉

+ 𝐶𝑜𝑒𝑓3 (𝑥) 𝐿

(

𝑘𝑚 , ℎ𝑏

0.25

0.50

0.75

1.00

X

(b)

′′ (𝜉) sin (𝑛𝜋𝜉∕𝐿)d𝜉 ∞ ∑ ∞ ∑

0.00

Fig. 2. The dimensionless displacement W (X; X0 , Y0 ) (a) and the dimensionless temperature distribution T (X, Y; X0 , Y0 ) (b) of the simple supported Euler– Bernoulli beam. Data are for Ω′ = 0.0, cs = ca = 0 and (X0 , Y0 ) = (0.5, 0.5).

) (31)

5. Analytical solutions of the coupled thermos-electro-elastic vibrations of the Euler–Bernoulli piezoelectric laminated beams

∞ ∑ ∞ ∑ ) ( ) 𝐿( 1 − (−1)𝑛 𝐹01 (𝑚, 𝑛)Φ𝑚 𝑘𝑚 , ℎ𝑏 𝑛𝜋 𝑚=1 𝑛=1

On the basis of the foregoing analysis in Section 5, the differential system Eq. (15) and the integral equation Eq. (29) are in fact equivalence. Both can be used to express the coupled thermoelectroelastic vibrations of laminated beams. By thus solving the Eq. (29) can fulfill our purpose. As a matter of fact, the solution of the integral equation Eq. (29) can be derived using the kernel function method [44].

× ∫0 𝑊 ′′ (𝜉) sin (𝑛𝜋𝜉∕𝐿)d𝜉 ∞ ∑ ∞ ∑ ( ) 𝐹01 (𝑚, 𝑛)𝑓𝐷1 (𝑛, 𝑥)Φ𝑚 𝑘𝑚 , ℎ𝑏 + 𝐶𝑇 2 𝑚=1 𝑛=1 𝐿

× ∫0 𝑊 ′′ (𝜉) sin (𝑛𝜋𝜉∕𝐿)d𝜉 + 𝑁𝑜𝑛(𝑥; 𝑥0 , 𝑦0 )

Table 1 The dimensionless displacement amplitude W (0.5 0.5; 0.5) and temperature T (0.5, 0.0; 0.5, 0.5) of the simply supported piezoelectric laminated beam. Orders of the present solutions

W(0.5; 0.5, 0.5) × 109

T(0.5, 0.0; 0.5, 0.5) × 102

1 3 5 7 10 20 30 40 FEM solution Diff3 (%)

4.1278 5.3444 5.5453 5.6246 5.6667 5.7183 5.7271 5.7298 5.6762 0.94%

1.1041 1.7916 2.6283 3.1885 3.6966 4.9114 5.2726 5.3907 5.3781 0.23%

Diff1 (%)

Diff2 (%)

22.764% 3.6228% 1.4098% 0.7429% 0.9024% 0.1537% 0.0471%

38.374% 31.834% 17.569% 13.745% 24.734% 6.8505% 2.1908%

Note: 𝐷𝑖𝑓 𝑓1 = |(𝑊 𝑗 − 𝑊 𝑖 )∕𝑊 𝑗 |, 𝐷𝑖𝑓 𝑓2 = |(𝑇 𝑗 − 𝑇 𝑖 )∕𝑇 𝑗 | and 𝐷𝑖𝑓 𝑓3 = |(∗𝐹 𝐸𝑀 − ∗40 )∕∗𝐹 𝐸𝑀 |. 360

X. Zhao, F.J.N. Iegaink and W.D. Zhu et al.

International Journal of Mechanical Sciences 156 (2019) 355–369

Fig. 3. The dimensionless voltage | V/(Y0 𝜔2 ) | as a function of the frequency (a) and the dimensionless displacement as a function of the frequency (b). Data are for R l = 106 Ω. Fig. 4. The dimensionless coupled displacement amplitude W (X; X0 , Y0 ) (a) and thermo-electrical combined coupling effect W1 (X; X0 , Y0 ) (b) for different heating positions of the concentrated heat load.

More specifically, taking derivatives to the Eq. (31) twice and multiplying sin (n 𝜋 𝜉/L) on the two sides of Eq. (31) and integrating once from zero to L leads to the following algebraic equations ( ) (32) 𝐈 − 𝐶𝑇 1 𝐀 − 𝐁1 − 𝐁2 − 𝐁3 − 𝐶𝑇 2 𝐂 𝐗 = 𝐛 where ( ) 𝐗𝑇 = 𝑋 1 , 𝑋 2 , ⋯ , 𝑋 𝑟2 , 𝑋 𝑛 =

𝐿

∫0

( 𝑊 ′′ (𝜉) sin

) 𝑛𝜋𝜉 d𝜉 𝐿

The specific terms in these matrices in the Eq. (34) are listed in Appendix B. In addition, the vector XT means the transposition of the vector X. Solving the algebraic equations Eq. (32), we can obtain the coupled parameters Xn . Then substituting the Xn into the Eq. (31) give rise to the explicit solution of the transversal displacement

(33)

and ⎛ 𝑎11 ⎜𝑎 𝐀 = ⎜ 21 ⎜⋯ ⎜𝑎 ⎝ 𝑟2 1 ⎛ 𝑏11 ⎜ 𝑏2 𝐁2 = ⎜ 21 ⎜⋯ ⎜ 2 ⎝𝑏𝑟 1 2

2

⎛ 𝑐11 ⎜𝑐 𝐶 = ⎜ 21 ⎜⋯ ⎜𝑐 ⎝ 𝑟2 1

𝑎12 𝑎22 ⋯ 𝑎𝑟2 2 𝑏212 𝑏222

⋯ 𝑏2𝑟 2 2

𝑐12 𝑐22 ⋯ 𝑐𝑟2 2

⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯ ⋯

⎛ 𝑏11 𝑎1𝑟2 ⎞ ⎜ 𝑏1 𝑎2𝑟2 ⎟⎟ , 𝐁1 = ⎜ 21 ⎜⋯ ⋯ ⎟ ⎜ 1 𝑎𝑟2 𝑟2 ⎟⎠ ⎝𝑏𝑟2 1

𝑏112

⋯

⋯ 𝑏1𝑟 2

⋯ ⋯

3 ⎛ 𝑏11 ⎜ 𝑏3 ⎜ 21

𝑏312 𝑏322

1

𝑏21𝑟 ⎞ 2 𝑏22𝑟 ⎟ 2 ⎟, 𝐁 = 3 ⎜⋯ ⋯ ⎟ ⎟ ⎜ 3 2 𝑏𝑟 𝑟 ⎠ ⎝𝑏𝑟2 1 2 2

𝑏122 2

⋯ 𝑏3𝑟 2 2

⋯

⋯ ⋯ ⋯ ⋯

𝑊 (𝑥; 𝑥0 , 𝑦0 ) = 𝑊1 (𝑥)+𝑊2 (𝑥; 𝑥0 , 𝑦0 ) ) ∞ ∑ ∞ ( 2 2 ∑ = 𝐶𝑇 1 − 𝑛𝐿𝜋2 𝐹01 (𝑚, 𝑛)𝑓𝐷1 (𝑛, 𝑥)𝐹11 (𝑚)𝑋𝑛

𝑏11𝑟 ⎞ 2 𝑏12𝑟 ⎟ 2 ⎟, ⋯ ⎟ ⎟ 𝑏1𝑟 𝑟 ⎠

𝑚=1 𝑛=1

+ 𝐶𝑜𝑒𝑓1 (𝑥)

2 2

𝑏31𝑟 2 𝑏32𝑟 2

⎞ ⎟ ⎟, ⋯ ⎟ ⎟ 𝑏3𝑟 𝑟 ⎠ 2 2

+ 𝐶𝑜𝑒𝑓2 (𝑥) (34)

+

𝑚=1 𝑛=1 ∞ ∑ ∞ ∑ 𝑚=1 𝑛=1 ∞ ∑ ∞ ∑

− 𝑛𝐿𝜋2

2 2

)

𝐹01 (𝑚, 𝑛)𝑓𝐷 (𝑛)𝐹11 (𝑚)𝑋𝑛

( ) 𝐹01 (𝑚, 𝑛)𝑓𝐷 (𝑛)Φ𝑚 𝑘𝑚 , ℎ𝑏 𝑋𝑛 (

)

( ) 𝑘𝑚 , ℎ𝑏 𝑋𝑛

𝐿 1 − (−1)𝑛 𝐹01 (𝑚, 𝑛)Φ𝑚 𝑛𝜋 𝑚=1 𝑛=1 ∞ ∑ ∞ ( ) ∑ 𝐶𝑇 2 𝐹01 (𝑚, 𝑛)𝑓𝐷1 (𝑛, 𝑥)Φ𝑚 𝑘𝑚 , ℎ𝑏 𝑋𝑛 𝑚=1 𝑛=1

+ 𝐶𝑜𝑒𝑓3 (𝑥)

𝑐1𝑟2 ⎞ ( ) 𝑐2𝑟2 ⎟⎟ 𝐛 = 𝑏1 , 𝑏2 , ⋯ , 𝑏𝑟2 . ⋯ ⎟ 𝑐𝑟2 𝑟2 ⎟⎠

∞ ∑ ∞ ( ∑

+ 𝑁𝑜𝑛(𝑥; 𝑥0 , 𝑦0 )

(35) 361

X. Zhao, F.J.N. Iegaink and W.D. Zhu et al.

International Journal of Mechanical Sciences 156 (2019) 355–369

Fig. 5. The dimensionless electro-displacement combined coupling effect T1 (X, Y; X0 , Y0 ) of the metallic layer for different heating positions of the concentrated heat load.

where W2 (x; x0 , y0 ) = Non (x; x0 , y0 ) and W1 (x; x0 , y0 ) is the double sum terms embrace the coupled parameters X n . It can be shown in the Eq. (35), the total displacement W (x; x0 , y0 ) embraces two terms W1 (x) and W2 (x; x0 , y0 ). More specifically, the term W1 (x) denotes the combined thermoelectric coupling of the total displacement W (x; x0 , y0 ) while the W2 (x; x0 , y0 ) stands for the uncoupled displacement solutions. Then according to the Eq. (28), the voltage expression can be also derived

where V1 is equal to the terms contain the coupled parameters Xn , and the other terms are V2 . Furthermore, based on the Eq. (18), it can arrive at the temperature field expression ( ) ( ) 𝑇 𝑥, 𝑦; 𝑥0 , 𝑦0 = 𝑇1 (𝑥, 𝑦) + 𝜆−1 𝑇 𝑞0 𝑇2 𝑥, 𝑦; 𝑥0 , 𝑦0 , 𝑇1 (𝑥, 𝑦) = (

[

∞ ∑ ∞ ∑

( 2

𝐹0 𝑚, 𝑛; 𝑥0 , 𝑦0

)

(37)

( ) sin (𝑛𝜋𝑥∕𝐿)Φ𝑚 𝑘𝑚 , 𝑦 .

It can be seen in the expressions of the voltage and the temperature that as well as the displacement, the voltage and the temperature also consist of two parts. In fact, the terms V1 and T1 (x, y), which have the coupled parameters X n , stands for combined coupling effect thermo-displacement of the total voltage and the combined coupling effect electro-displacement of the total temperature, respectively. In addition, V2 (x0 , y0 ) and T2 (x, y; x0 , y0 ) are the uncoupled terms.

]

1−𝐶𝐷3 ⋅𝐶𝐷2 ⋅𝐶𝐺1

∑( ∑( ) ) 𝜆−1 𝑚, 𝑛; 𝑥0 , 𝑦0 + 𝐶𝑇 2 𝐶𝐷3 𝑚, 𝑛; 𝑥0 , 𝑦0 𝑇 𝑞0 𝐶𝑇 1 𝐶𝐷3 1 3 ] ∑( ) +𝑝𝑇 𝑏𝑝 𝑚, 𝑛; 𝑥0 , 𝑦0

+

)

( ) 𝐹01 (𝑚, 𝑛) sin (𝑛𝜋𝑥∕𝐿)Φ𝑚 𝑘𝑚 , 𝑦 𝑋𝑛 ,

𝑚=1 𝑛=1

∞ ∑ ∞ ∑ ) ( ) 𝐿( +𝑝𝑇 𝑏𝑝 1 − (−1)𝑛 𝐹01 (𝑚, 𝑛)Φ𝑚 𝑘𝑚 , ℎ𝑏 𝑋𝑛 𝑛𝜋 𝑚=1 𝑛=1

𝑉 = 𝑉1 + 𝑉2 = [

𝑚=1 𝑛=1

𝑇2 𝑥, 𝑦; 𝑥0 , 𝑦0 =

( 2 2) 𝑛 𝜋 𝐶𝑇 1 𝐶𝐷3 − 𝐹01 (𝑚, 𝑛)𝑓𝐷 (𝑛)𝐹11 (𝑚)𝑋𝑛 𝐿2 𝑚=1 𝑛=1 ∞ ∑ ∞ ∑ ( ) +𝐶𝑇 2 𝐶𝐷3 𝐹01 (𝑚, 𝑛)𝑓𝐷 (𝑛)Φ𝑚 𝑘𝑚 , ℎ𝑏 𝑋𝑛 ∞ ∑ ∞ ∑

𝑚=1 𝑛=1

∞ ∑ ∞ ∑

6. Numerical results and discussions To this point, the majority of this review has concentrated on ways to obtain the analytical solutions of our problem. It is therefore important to realize some numerical results. We will start by considering a simple supported piezoelectric laminated beam with a PZT 5A &5H layer and

5

1 − 𝐶𝐷3 ⋅ 𝐶𝐷2 ⋅ 𝐶𝐺1 (36) 362

X. Zhao, F.J.N. Iegaink and W.D. Zhu et al.

International Journal of Mechanical Sciences 156 (2019) 355–369

static deflection at the middle section x0 = L/2 of the EB subject to a unit force at x0 = L/2. The T0 is the initial temperature. The 𝜁 1 r and 𝜁 2 r respectively denote the air damping and material damping. The 𝛽 is the height-to-length ratio of the beam, and the 𝛽 1 is the composing proportion of the PZT layer and the aluminum layer. The material parameters used in the numerical example are listed in the following. It is assumed that all the material properties are constant over the temperature range studied in this work. Piezoelectric: ( ) ( ) 𝐸𝑝 = 6.7 × 1010 N∕m2 , 𝜌 = 7800 kg∕m3 , 𝑑31 = −1.90 × 10−10 (m∕V), 𝜀𝑠33 = 15.93 × 10−9 (F∕m), 𝑝𝑇 = −2.5 × 10−5 (CK −1 m−2 ); Aluminum: 𝐸 = 7.2 × 1010 (N∕m2 ), 𝜌 = 2700(kg∕m3 ), 𝑐𝑝 = 880(J∕kg ⋅ K), 𝜆𝑇 = 237(𝑊 ∕(m ⋅ K )), 𝛼𝑇 = 2.30 × 10−5 (1∕K). In this section, the following specific thermal boundary condition will be considered in the following part ( ) d𝐺𝑇 𝑥, ℎ𝑎 d𝑧

( ) d𝐺𝑇 𝑥, ℎ𝑏 ( ) + 𝜇1 𝐺𝑇 𝑥, ℎ𝑎 = 𝜇1 𝑇0 , = 0, d𝑧 𝐺𝑇 (0, 𝑧) = 𝐺𝑇 (𝐿, 𝑧) = 0.

(39)

In the Eq. (39), the relative heat transfer coefficient 𝜇1 = 𝛼 1 /𝜆 T , 𝛼 1 = 28,000 (W/(m2 ∙K)). Moreover, the shear modified coefficient is 𝜅 = 0.83, the heat flux density Q0 = 50,000 (W/m3 ), the electric resistance R l = 100 (Ω), 𝛽 = 0.1 and 𝛽 1 = 0.25. For simplicity, except the last numerical example, the damping effects are not considered in the other examples. Besides the first and the fourth examples, the other examples take the external heat load frequency Ω′ = 0.5. Since the external load frequency is identified, the two types of damping can be simply expressed by 𝜁 1 and 𝜁 2 .

6.1. Verification of the present solution In this subsection, the validation of the present solution will be verified by FEM results, and the convergence of the present series solutions will be also shown. Since the dynamic calculation could consume a tremendous of time, we simplicity the work by letting the frequency of the external heat load equal to zero i.e. Ω′ = 0, and an isotropic uniform metallic beam, which has no piezoelectric layer, is considered in this example to simplify the verification process. Based on these degenerations, the coupled thermos-electro-elastic vibration problem degenerates to an uncoupled static thermal stress problem. Steel is taken as the material in this validation example especially to compare with the numerical results in Ref. [35]. As a first step, the convergence of the present solutions will be discussed. For simplicity, we let the order r1 equal to r2 . As shown in Table 1, as the order of the double series r1 (r2 ) increasing, the differences between the displacement amplitudes W (0.5; 0.5, 0.5) are decreasing. Meanwhile, for the temperature field, the temperature value T (0.5, 0.0; 0.5, 0.5) is also decreasing with the order r1 (r2 ). These facts can be used to validate the convergence of current solutions. Besides, as it has been found in Ref. [35], in Euler–Bernoulli case, it can also find that the convergent speed of the displacement is much faster than the temperature. To be more specific, it can be shown in Table 1 that for the displacement, when the order r1 (r2 ) = 5, the relative difference Diff1 = 3.6228% which is smaller than 5%. However, in the temperature case, the relative difference Diff2 is smaller than 5% until the order r1 (r2 ) = 40. Physically speaking, the velocity of the heat propagation is much slower than that of the elastic wave. Comparing with the results for the Timoshenko case reported in Ref. [35], the convergence of the temperature is more difficult. More specifically, for

Fig. 6. The dimensionless coupled displacement amplitude W (X; X0 , Y0 ) (a) and the thermo-electrical combined coupling effect W1 (X; X0 , Y0 ) (b) of a laminated beam for various electric resistances R l .

an aluminum layer. A concentrated heat load, whose heat flux density is Q0 , is exerted to the laminated beam at position (x0 , y0 ) i.e. q (x, y, t) = Q0 𝛿 (x – x0 ) 𝛿 (y – y0 ) ei𝜔 t . We introduce the following dimensionless quantities are introduced in this section. ( ) 𝑊 (𝑥; 𝑥0 , 𝑦0 ) 𝑦 𝑥 𝜔 𝑉 ,Y= , Ω′ = , 𝑊 X;X0 ,Y0 = , 𝑉𝐷 = , 𝐿 ℎ𝑠 𝜔0 𝐿 𝐶𝐷3 (Ω0 ) ( ) ( ) 𝑇 𝑥, 𝑦; 𝑥0 , 𝑦0 𝑉1 𝑉𝐷1 = , 𝑇 X, Y;X0 ,Y0 = , 𝐶𝐷3 (Ω0 ) 𝑇0 ( ) 𝑊 (𝑥; 𝑥 , 𝑦 ) 𝑐 𝐼𝜔 𝑐 𝜁𝑟 = 𝜁1𝑟 + 𝜁2𝑟 = 𝑠 𝑟 + 𝑎 , 𝑊1 X;X0 ,Y0 = 1 𝑠 0 0 , 2𝐸𝐼 2𝑚𝜔𝑟 𝑊max X=

𝑇1 (X, Y) =

ℎ𝑝 𝑇1 (𝑥, 𝑦) ℎ , 𝛽 = , 𝛽1 = 𝑇0 𝐿 ℎ𝑠

(38)

where the symbol 𝜔0 = 𝜋 2 (EI/𝜌A)0.5 /L2 denotes the first order natural frequency of the Euler–Bernoulli beam. The Ws max = L/(48EI) is the 363

X. Zhao, F.J.N. Iegaink and W.D. Zhu et al.

International Journal of Mechanical Sciences 156 (2019) 355–369

Fig. 7. The dimensionless electro-displacement combined coupling effect T1 (X, Y; X0 , Y0 ) of the metallic layer for various electric resistances R l .

the Timoshenko beam, the relative difference Diff2 = 3.1980% < 5% just needs the order r1 (r2 ) = 10. However, for the Euler–Bernoulli beam, the same purpose Diff2 < 5% is satisfied when the order r1 (r2 ) = 40. The FEM solutions are utilized to validate the present solutions. It can be revealed in Table 1 that the differences between the present solutions and the FEM solutions are smaller than 1%, whether it is for the displacement or for the temperature. This fact can be used to verify the present results. Furthermore, in Fig. 2(a), it is shown that the FEM displacement is of agreement with the present displacement, which can be also applied to verify the present solutions. Moreover, the global temperature distribution is also furnished in Fig. 2(b). Another degeneration of the coupled thermo-electro-elastic vibration equation of the laminated beam can be used to verify the voltage solution. In this case, the heat transfer process is intersected and the heat transfer coefficient 𝛼 T = 0 and the pyroelectric parameter pT = 0. From the Eq. (15), it can be seen that the parameter CT 1 = CT 2 = 0 and the CD 4 = 0 which delete the thermal coupled terms of the system. In view of utilizing Erturk and Inman’ results [12], a harmonic force F (x, t) = F (x) ei𝜔 t is exerted to the laminated beam. Thus the Eq. (15) becomes 𝑊 ′ + 𝐶𝐷1 𝑊 + 𝐶𝐷2 ⋅ 𝑉 [𝛿 ′ (𝑥 − 𝑥1 ) − 𝛿 ′ (𝑥 − 𝑥2 )] = 𝐹 (𝑥),

Δ𝑇 + iΩ𝐶𝑇 3 𝑇 = −𝜆−1 𝑇 𝑞 (𝑥, 𝑦) 𝑉 = 𝐶𝐷3

𝑥2

∫𝑥1

𝑊 ′′ 𝑑𝑥

(40b) (40c)

Since the thermal coupled terms vanish, the heat transfer process Eq. (40b) cannot induce thermal stress to influence the deformation of the beam. Therefore, the Eq. (40b) can be ignored. Then the Eq. (40) can be written as 𝑊 ′ + 𝐶𝐷1 𝑊 + 𝐶𝐷2 ⋅ 𝑉 [𝛿 ′ (𝑥 − 𝑥1 ) − 𝛿 ′ (𝑥 − 𝑥2 )] = 𝐹 (𝑥), 𝑉 = 𝐶𝐷3

𝑥2

∫𝑥1

𝑊 ′′ 𝑑𝑥,

(41a) (41b)

which is just the steady-state vibration form in the work of Erturk and Inman [12]. Employing the Green’s functions in Section 4, we can get the closed-form solution of the voltage 𝑥

𝑉 =

𝑥

364

1

1 + 𝐶𝐷2 ∫𝑥 2 [𝐺𝐷2 ′′ (𝑥; 𝑥1 ) − 𝐺𝐷2 ′′ (𝑥; 𝑥2 )]d𝑥 1

(40a)

𝐿

𝐶𝐷3 ∫𝑥 2 ∫0 𝐹 (𝜉)𝐺𝐷1 ′′ (𝑥; 𝜉)𝑑𝜉d𝑥

(42)

X. Zhao, F.J.N. Iegaink and W.D. Zhu et al.

International Journal of Mechanical Sciences 156 (2019) 355–369

Fig. 9. The dimensionless coupled displacement amplitude W (X; X0 , Y0 ) (a) and the thermo-electrical combined coupling effect W1 (X; X0 , Y0 ) of the laminated beam for varying air damping 𝜁 1 .

Fig. 8. The dimensionless coupled voltage | VD | of a simple supported piezoelectric laminated beam (a) and a cantilever piezoelectric laminated beam | VD | (b) for different convective heat transfer coefficient 𝛼 1 .

amplitudes W2 (X; X0 , Y0 ) is not illustrated in the Fig. 4. The reason is that W2 (X; X0 , Y0 ) actually equal to W (X; X0 , Y0 )−W1 (X; X0 , Y0 ). Thus it is easy to estimate the value of W2 (X; X0 , Y0 ) through W (X; X0 , Y0 ) and W1 (X; X0 , Y0 ). It can be shown in the Fig. 4(a) that if the heating position is different, the coupled displacements have a negative or positive value and the relations between the coupled and uncoupled displacements are also various. In more specific terms, when the (x0 , y0 ) = (0.5L, ha /2), the W (X; X0 , Y0 ) < 0, while when the (x0 , y0 ) = (0.5L, 0.0) and (x0 , y0 ) = (0.5L, ha ), the W (X; X0 , Y0 ) > 0. The value of W (X; 0.5L, hb ) has both positive and negative part. For the thermo-electrical combined coupling effect W1 (X; X0 , Y0 ), from the Fig. 4(b), it is seen that when the (X0 , Y0 ) = (0.5, ha /2) and (X0 , Y0 ) = (0.5, hb ), the W1 (X; X0 , Y0 ) is convex function. However, the W1 (X; X0 , Y0 ) becomes concave function as the (X0 , Y0 ) = (0.5, ha ) and (X0 , Y0 ) = (0.5, 0.0). Physically, the thermo-electrical combined coupling effect W1 (X; X0 , Y0 ) is triggered by the temperature field T1 (X, Y; X0 , Y0 ) which means the electro-displacement combined coupling effect. The Fig. 5 displays the dimensionless electro-displacement combined coupling effect T1 (X,

Using the cantilever boundary conditions and the harmonic based exciting F (x) = Y0 𝜇𝜔2 where the Y0 denotes the amplitude of the based exciting, we can obtain the value of the voltage. In this example, the damping effects are not considered. It is shown in Fig. 3 that the present voltages and displacements agree with the convergent solutions of the Ref. [12] perfectly. These facts provide verifications of the present solutions on another side. 6.2. Influence of the heating positions to the displacement and temperature solutions In this subsection, the influence of the heating positions of a concentrated heat load q (x, y, t) = Q0 𝛿 (x – x0 ) 𝛿 (y – y0 ) ei𝜔 t to the present solutions will be discussed. The Fig. 4 shows the dimensionless coupled displacement amplitude W (X; X0 , Y0 ) and thermo-electrical combined coupling effect W1 (X; X0 , Y0 ) under different heating positions. The uncoupled displacement 365

X. Zhao, F.J.N. Iegaink and W.D. Zhu et al.

International Journal of Mechanical Sciences 156 (2019) 355–369

Fig. 10. The dimensionless electro-displacement combined coupling effect T1 (X, Y; X0 , Y0 ) of the metallic layer for various air damping 𝜁 1 .

Y; X0 , Y0 ) of the metallic layer under different heating positions. Since the electrical effect is considered, the temperature field is more complex than the thermoelastic vibration case [35]. More specifically, when the (X0 , Y0 ) = (0.5, ha ) and (X0 , Y0 ) = (0.5, 0.0), the cold-contraction effects of the “ice core regions” on the upper and lower surfaces can be offset each other. Therefore, the W1 (X; X0 , Y0 ) is actually identified by the “hot core region” located under the neutral surface. According to the location of the “hot core region” and the heat-expansion of it, it can be known the curves of these two typical cases are concave. The same analysis can be also applied to the situations (X0 , Y0 ) = (0.5, ha /2) and (X0 , Y0 ) = (0.5, hb ). For the (X0 , Y0 ) = (0.5, ha /2) and (X0 , Y0 ) = (0.5, hb ), the determiner changes to the “ice core region” under the neutral surface. Thus the shapes of the curves are convex. The talked “ice or hot core region” phenomena have been used in Refs [34,35] to explain some problems.

placement W (X; X0 , Y0 ) and the thermo-electrical combined coupling effect W1 (X; X0 , Y0 ) for various electric resistances R l . From the Fig. 6(a), it is seen that when the electric resistances are smaller Rl = 102 Ω (short circuit [12]), the coupled displacement W (X; X0 , Y0 ) > 0. However, as the electric resistance increases, the W (X; X0 , Y0 ) becomes negative for the cases R l = 103 Ω, R l = 104 Ω and Rl = 105 Ω. For the thermo-electrical combined coupling effect W1 (X; X0 , Y0 ) in the Fig. 6(b), globally, the amplitude of W1 (X; X0 , Y0 ) decreases with the resistances R l . More specifically, for the R l = 102 Ω and R 3 4 l = 10 Ω, the W1 (X; X0 , Y0 ) < 0. Nonetheless, for the Rl = 10 Ω and 5 Rl = 10 Ω, the W1 (X; X0 , Y0 ) has both positive and negative parts, and the positive value is becoming larger as the electric resistance increases. As a matter of fact, the phenomena of the W1 (X; X0 , Y0 ) can be explained by the coupling effect temperature T1 (X, Y; X0 , Y0 ). Fig. 7 plots the dimensionless electro-displacement combined coupling effect T1 (X, Y; X0 , Y0 ) of the metallic layer for various electric resistances Rl . From the physical point of view, it can be revealed in Fig. 7 that for the Rl = 102 Ω, the shape of the curve is determined by the “hot cores” located under the neutral surface. The heat-expansion effect concaves the W1 (X; X0 , Y0 ). However, as the electric resistance becoming larger, the T1 (X, Y; X0 , Y0 ) exhibits different forms for the Rl = 103 Ω, Rl = 104 Ω and Rl = 105 Ω. To be more specific terms, the “ice core” under the neutral surface make the middle area of the beam become convex gen-

6.3. Influence of the electric resistances to the displacement and temperature solutions The essential circuit parameter—electric resistance will be discussed in this subsection. It strives to reveal the effect of the electric resistances on the temperature and displacement of the beam. In this numerical example, the heating position is set to (X0 , Y0 ) = (0.5, 0.0). The Fig. 6 exhibits the dimensionless coupled dis366

X. Zhao, F.J.N. Iegaink and W.D. Zhu et al.

International Journal of Mechanical Sciences 156 (2019) 355–369

erally. Moreover, as the temperature increases around the “ice core”, two “hot cores” appear close to the middle “ice core”, which make the perimeter areas concave generally. For engineering application, the specific temperature distribution under the open circuit or short circuit condition could be used as a benchmark to check the circuit condition. 6.4. Influence of the convective heat transfer coefficient to the voltage solutions This subsection will discuss the effect of the heat parameter– convective heat transfer coefficient on the voltage solutions. The corresponding discussion is related to the energy harvesting problems of the system. As same as Section 6.3, the heating position is also (X0 , Y0 ) = (0.5, 0.0). The Fig. 8 shows the dimensionless coupled voltage | VD | of a simply supported piezoelectric laminated beam and a cantilever piezoelectric laminated beam | VD | for different convective heat transfer coefficient 𝛼 1 . In present studies for the piezoelectric energy harvester, most of them are the cantilever piezoelectric laminated beams. Therefore, two different boundary conditions are discussed in this subsection. Initially, for the simply supported one, it can be displayed from the Fig. 8(a) that the maximum value of the dimensionless coupled voltage | VD | = 46.77 at 𝛼 1 = 2.0 × 105 . To be more specific, when the (𝛼 1 , Ω′) = (2.0 × 104 , 4.92), the coupled voltage reaches the maximum value | VD | = 46.77, and the best quality bandwidth is also aroused at the 𝛼 1 = 2.0 × 104 . Considering the bandwidth and the extremum value of the | VD |, the best energy harvesting coefficient should take the 𝛼 1 = 2.0 × 104 . Furthermore, for the cantilever case, from the Fig. 8(b), it can be also found that the 𝛼 1 = 2.0 × 104 is definitely the best energy harvesting parameter. At the point (𝛼 1 , Ω′) = (2.0 × 104 , 6.00), the peak value of the dimensionless coupled voltage | VD | = 211.47, which is much larger than the simply supported case. This is proving that cantilever beam is much better than a simply supported beam to be piezoelectric energy harvesters. Incidentally, at the point (𝛼 1 , Ω′) = (1.5 × 104 , 6.00), the voltage | VD | = 190.56. Therefore, the parameter 𝛼 1 = 1.5 × 104 is also a good choice. 6.5. Influence of the two damping parameters to the present solutions In this subsection, the effects of the air damping and the material damping on the present solutions will be discussed. Varying from the above subsection, we take the heating position (X0 , Y0 ) = (0.5, ha ) in the following numerical examples. The Fig. 9 plots the dimensionless coupled displacement W (X; X0 , Y0 ) (a) and the thermo-electrical combined coupling effect W1 (X; X0 , Y0 ) for varying air damping coefficient 𝜁 1 . As expected, the coupled displacement amplitudes W (X; X0 , Y0 ) and W1 (X; X0 , Y0 ) are decreasing with the air damping coefficient 𝜁 1 . As a matter of fact, the displacement W1 (X; X0 , Y0 ) is induced by the electro-displacement combined coupling effect T1 (X, Y; X0 , Y0 ). Fig. 10 displays the dimensionless electro-displacement combined coupling effect T1 (X, Y; X0 , Y0 ) of the metallic layer for various air damping 𝜁 1 . It can be revealed that a “hot core” and an “ice core” are formed in Fig. 10. Since the heat-expansion of the “hot core” under the neutral surface concaves the W1 (X; X0 , Y0 ) but the cold-contraction effect of the “ice core” on the lower surface become convex, the shape and the amplitude of the curve W1 (X; X0 , Y0 ) is determined by the competition of the two opposite effects. As the air damping 𝜁 1 increasing, the middle area of the “hot core” become concave and move down toward the bottom. The temperature of the “hot core” is reducing generally. On the contrary, the temperature of the “ice core” become colder and colder. These facts show that the W1 (X; X0 , Y0 ) is concave firstly and then become convex generally in the middle region of the curve. The amplitude is also reducing in this process.

Fig. 11. The dimensionless coupled displacement amplitude W (X; X0 , Y0 ) (a) and the thermo-electrical combined coupling effect W1 (X; X0 , Y0 ) (b) of the laminated beam for different material damping 𝜁 2 .

For the material damping effect, Fig. 11 shows the dimensionless coupled displacement amplitude W (X; X0 , Y0 ) and the thermo-electrical combined coupling effect W1 (X; X0 , Y0 ) of the laminated beam for different material damping 𝜁 2 . It can be seen in the Fig. 11(a) that the coupled displacement is decreasing with the material damping observably. From Fig. 11(b), it can be obtained the amplitude of the W1 (X; X0 , Y0 ) increase with the material damping 𝜁 2 . More specifically, the enhancement of the W1 (X; X0 , Y0 ) is also triggered by the electro-displacement combined coupling effect T1 (X, Y; X0 , Y0 ). Fig. 12 is the dimensionless electro-displacement combined coupling effect T1 (X, Y; X0 , Y0 ) of the metallic layer for varying material damping 𝜁 2 . It can be known that the shape of the curve W1 (X; X0 , Y0 ) is determined by the competition of the “hot core” under the neutral surface and the “ice core” on the lower surface. As the same with the previous analysis, when the material damping increasing, the middle region of the “hot core” become convex and move up towards the neutral surface generally, and the temperature value of the “hot core” increases with the material damping. These phenomena reveal the changing of the W1 (X; X0 , Y0 ). 367

X. Zhao, F.J.N. Iegaink and W.D. Zhu et al.

International Journal of Mechanical Sciences 156 (2019) 355–369

Fig. 12. The dimensionless electro-displacement combined coupling effect T1 (X, Y; X0 , Y0 ) of the metallic layer for varying material damping 𝜁 2 .

7. Conclusions

Acknowledgments

In this paper, the coupled thermos-electro-elastic vibrations of the two-layer laminated Euler–Bernoulli beams with a simply supported boundary condition are investigated. The classical electric circuit model is developed based on the piezothermoelastic constitutive relations. Green’s function method is utilized to transform the coupled differential equations to an uncoupled integral equation. The explicit solutions of the displacement, temperature and voltage are derived. On the basis of these analytical expressions, some numerical examples have been conveniently done to analyze the influences of some essential parameters. We then used the FEM solutions to verify the current results. Another verification method is to compare the present solutions with others results from some references. Through analyzing the effect of the electric resistances to the temperature field T1 (X, Y; X0 , Y0 ) and the displacement W1 (X; X0 , Y0 ), it can be concluded that the specific temperature distribution under the open circuit condition could be used as a benchmark to check the circuit condition. Furthermore, according to the influence of the heat convection coefficient to the voltage, the optimized heat transfer coefficients are identified. These results may give some assistance to analyze and designing smart sensors or actuators and solving some engineering thermal problems.

The work was supported by the National Natural Science Foundation of China (nos. 11702230, 11772100, and 11372257) and Program for the risk assessments and control of the coastal pipe and its terrestrial terminal (no. 2016YFC0802305) and for Sichuan youth scientific and technological innovation research team of engineering structural safety assessment and disaster prevention technology (2019JDTD0017). The supports from Chong Qing Municipal Solid Waste Resource Utilization & Treatment Collaborative Innovation Center (grant number Shljzyh2017-007) and Southwest petroleum university bridge safety assessment youth science and technology innovation team (2018CXTD07). Appendix A The expressions of the Σi (i = 1, 2 … 10) in the text are listed here ) ∞ ∑ ∞ ( ∑ ∑ ( ) 𝑛2 𝜋 2 (𝑚, 𝑛; 𝑥0 , 𝑦0 ) = − 𝐹02 𝑚, 𝑛; 𝑥0 , 𝑦0 𝑓𝐷 (𝑛)𝐹11 (𝑚), 𝐿2 1 𝑚=1 𝑛=1 ) ∞ ∑ ∞ ( ∑ ∑ 𝑛2 𝜋 2 (𝑚, 𝑛; 𝜉) = − 𝐹01 (𝑚, 𝑛)𝑓𝐷 (𝑛)𝐹11 (𝑚) sin (𝑛𝜋𝜉∕𝐿), 𝐿2 2 𝑚=1 𝑛=1

368

X. Zhao, F.J.N. Iegaink and W.D. Zhu et al.

∑

(𝑚, 𝑛; 𝑥0 , 𝑦0 ) =

𝑚=1 𝑛=1

3

∑

(𝑚, 𝑛; 𝜉) =

4

∑

∞ ∑ ∞ ∑

International Journal of Mechanical Sciences 156 (2019) 355–369

( ) ( ) 𝐹02 𝑚, 𝑛; 𝑥0 , 𝑦0 𝑓𝐷 (𝑛)Φ𝑚 𝑘𝑚 , ℎ𝑏 ,

∞ ∑ ∞ ∑

𝐹01 (𝑚, 𝑛)𝑓𝐷 (𝑛)Φ𝑚 𝑚=1 𝑛=1 ∞ ∑ ∞ ∑

(𝑚, 𝑛; 𝑥0 , 𝑦0 ) =

𝑚=1 𝑛=1

5

[8] Wang D, Liu J, Zhou D, Huang S. Using PVDF piezoelectric film sensors for in situ measurement of stayed-cable tension of cable-stayed bridges. Smart Mater Struct 1999;8:554. [9] Zi Y, Lin L, Wang J, Wang S, Chen J, Fan X, Yang PK, Yi F, Wang ZL. Triboelectric-pyroelectric-piezoelectric hybrid cell for high-efficiency energy-harvesting and self-powered sensing. Adv Mater 2015;27:2340–7. [10] Athenstaedt H, Claussen H, Schaper D. Epidermis of human skin: pyroelectric and piezoelectric sensor layer. Science 1982;216:1018. [11] Newnham RE, Skinner DP, Cross LE. Connectivity and piezoelectric–pyroelectric composite. Mater Res Bull 1978;13:525–36. [12] Erturk A, Inman DJ. A distributed parameter electromechanical model for cantilevered piezoelectric energy harvesters. J Vib Acoust 2008;130:1257–61. [13] Danesh-Yazdi AH, Elvin N, Andreopoulos Y. Green‫׳‬s function method for piezoelectric energy harvesting beams. J Sound Vib 2014;333:3092–108. [14] Dietl JM, Wickenheiser AM, Garcia E. A Timoshenko beam model for cantilevered piezoelectric energy harvesters. Smart Mater Struct 2010;19:055018. [15] Zhao X, Yang E, Li Y, Crossley W. Closed-form solutions for forced vibrations of piezoelectric energy harvesters by means of Green’s functions. J Intell Mater Syst Struct 2017;28:2372–87. [16] Guan M, Liao W-H. Design and analysis of a piezoelectric energy harvester for rotational motion system. Energy Convers Manage 2016;111:239–44. [17] Fan K, Chang J, Chao F, Pedrycz W. Design and development of a multipurpose piezoelectric energy harvester. Energy Convers Manage 2015;96:430–9. [18] Zhang J, Fang Z, Shu C, Zhang J, Zhang Q, Li C. A rotational piezoelectric energy harvester for efficient wind energy harvesting. Sens Actuators A 2017;262:123–9. [19] Xiong H, Wang L. Piezoelectric energy harvester for public roadway: on-site installation and evaluation. Appl Energy 2016;174:101–7. [20] Mindlin RD. Equations of high frequency vibrations of thermopiezoelectric crystal plates. Int J Solids Struct 1973;10:625–37. [21] Tzou HS, Bao Y. A theory on anisotropic piezothermoelastic shell laminates with sensor/actuator applications. J Sound Vib 1995;184:453–73. [22] Jiang JP, Li DX. A new finite element model for piezothermoelastic composite beam. J Sound Vib 2007;306:849–64. [23] Li P, Li C, Cong B. Piezoelectric-thermo-elastic coupling effect analysis for piezoelectric vibration energy harvester. Microsyst Technol 2018;24(9):3823–32. [24] Wang YQ, Zu JW. Nonlinear steady-state responses of longitudinally traveling functionally graded material plates in contact with liquid. Compos Struct 2017;164:130–44. [25] Wang YQ, Yang Z. Nonlinear vibrations of moving functionally graded plates containing porosities and contacting with liquid: internal resonance. Nonlinear Dyn 2017;90:1–20. [26] Du C, Li Y, Jin X. Nonlinear forced vibration of functionally graded cylindrical thin shells. Thin-Walled Struct 2014;78:26–36. [27] Komijani M, Gracie R. Nonlinear thermo-electro-mechanical dynamic behaviour of FGPM beams. Compos Struct 2016;150:208–18. [28] Larkin K, Abdelkefi A. Neutral axis modeling and effectiveness of functionally graded piezoelectric energy harvesters. Compos Struct 2019;213:25–36. [29] Amini Y, Emdad H, Farid M. Finite element modeling of functionally graded piezoelectric harvesters. Compos Struct 2015;129:165–76. [30] Cao Y, Huang H, Zhu ZH, Su S. Optimized energy harvesting through piezoelectric functionally graded cantilever beams. Smart Mater Struct 2019;28:025038. [31] Ebrahimi F, Salari E. Analytical modeling of dynamic behavior of piezo-thermo– electrically affected sigmoid and power-law graded nanoscale beams. Appl Phys A 2016;122:793. [32] Jung HJ, Song Y, Hong SK, Yang CH, Hwang SJ, Jeong SY, Sung TH. Design and optimization of piezoelectric impact-based micro wind energy harvester for wireless sensor network. Sens Actuators A 2015;222:314–21. [33] Kuo C-L, Lin S-C, Wu W-J. Fabrication and performance evaluation of a metal-based bimorph piezoelectric MEMS generator for vibration energy harvesting. Smart Mater Struct 2016;25:105016. [34] Zhao X, Zhao YR, Gao XZ, Li XY, Li YH. Green‫׳‬s functions for the forced vibrations of cracked Euler–Bernoulli beams. Mech Syst Signal Process 2016:155–75 s 68–69. [35] Zhao X, Yang EC, Li YH. Analytical solutions for the coupled thermoelastic vibrations of Timoshenko beams by means of green’s functions. Int J Mech Sci 2015;100:50–67. [36] Grover D, Sharma JN. Transverse vibrations in piezothermoelastic beam resonators. J Intell Mater Syst Struct 2012;23:77–84. [37] Guo XX, Wang ZM, Wang Y, Zhou YF. Analysis of the coupled thermoelastic vibration for axially moving beam. J Sound Vib 2009;325:597–608. [38] Manoach E, Ribeiro P, Coupled thermoelastic. large amplitude vibrations of Timoshenko beams. Int J Mech Sci 2004;46:1589–606. [39] Massalas CV, Kalpakidis VK. Coupled thermoelastic vibrations of a Timoshenko beam. Int J Eng Sci 1984;22:459–65. [40] Zhao X, Hu QJ, Crossley W, Du CC, Li YH. Analytical solutions for the coupled thermoelastic vibrations of the cracked Euler-Bernoulli beams by means of Green’s functions. Int J Mech Sci 2017;128:37–53. [41] Asmar NH. Partial differential equations with Fourier series and boundary value problems. Prentice Hall; 2012. [42] Abu-Hilal M. Forced vibration of Euler–Bernoulli beams by means of dynamic Green functions. J Sound Vib 2003;267:191–207. [43] Li XY, Zhao X, Li YH. Green’s functions of the forced vibration of Timoshenko beams with damping effect. J Sound Vib 2014;333:1781–95. [44] Corduneanu C. Integral equations and applications. Cambridge University Press; 1991.

( ) 𝑘𝑚 , ℎ𝑏 sin (𝑛𝜋𝜉∕𝐿),

) ( ) ( ) 𝐿( 1 − (−1)𝑛 𝐹02 𝑚, 𝑛; 𝑥0 , 𝑦0 Φ𝑚 𝑘𝑚 , ℎ𝑏 . 𝑛𝜋

∑

∞ ∑ ∞ ∑ ) ( ) 𝐿( (𝑚, 𝑛; 𝜉) = 1 − (−1)𝑛 𝐹01 (𝑚, 𝑛)Φ𝑚 𝑘𝑚 , ℎ𝑏 sin (𝑛𝜋𝜉∕𝐿), 𝑛𝜋 6 𝑚=1 𝑛=1 ) ∞ ∑ ∞ ( ∑ ∑ ( ) 𝑛2 𝜋 2 (𝑚, 𝑛; 𝑥0 , 𝑦0 , 𝑥) = − 𝐹02 𝑚, 𝑛; 𝑥0 , 𝑦0 𝑓𝐷1 (𝑛, 𝑥)𝐹11 (𝑚), 2 𝐿 7 𝑚=1 𝑛=1 ) ∞ ∑ ∞ ( ∑ ∑ 𝑛2 𝜋 2 (𝑚, 𝑛; 𝑥, 𝜉) = − 𝐹01 (𝑚, 𝑛)𝑓𝐷1 (𝑛, 𝑥)𝐹11 (𝑚) sin (𝑛𝜋𝜉∕𝐿), 𝐿2 8 𝑚=1 𝑛=1

∑

𝑚=1 𝑛=1

9

∑

∞ ∑ ∞ ∑

(𝑚, 𝑛; 𝑥0 , 𝑦0 , 𝑥) = (𝑚, 𝑛; 𝑥, 𝜉) =

∞ ∑ ∞ ∑ 𝑚=1 𝑛=1

10

( ) ( ) 𝐹02 𝑚, 𝑛; 𝑥0 , 𝑦0 𝑓𝐷1 (𝑛, 𝑥)Φ𝑚 𝑘𝑚 , ℎ𝑏 ,

( ) 𝐹01 (𝑚, 𝑛)𝑓𝐷1 (𝑥, 𝑛)Φ𝑚 𝑘𝑚 , ℎ𝑏 sin (𝑛𝜋𝜉∕𝐿). (A1)

In the Eq. (A1), ( ) ℎ 𝐹11 (𝑚) = ∫ℎ 𝑏 Φ𝑚 𝑘𝑚 , 𝑦 𝑦d𝑦, 𝑎

𝑓𝐷 (𝑛) =

𝐿 𝐿 ∫0 ∫0

𝐺𝐷1 ′′ (𝑥; 𝜉) sin (𝑛𝜋𝜉∕𝐿)d𝜉d𝑥,

(A2)

𝐿

𝑓𝐷1 (𝑛, 𝑥) = ∫0 𝐺𝐷1 (𝑥; 𝜉) sin (𝑛𝜋𝜉∕𝐿)𝑑𝜉. Appendix B The expressions of the matrix terms in Eq. (34) are listed in the following 𝑎𝑖𝑗 = 𝑏1𝑖𝑗 = 𝑏2𝑖𝑗 = 𝑏3𝑖𝑗

=

𝑐𝑖𝑗 = 𝑏𝑛 =

𝑟1 ( ∑ 𝑚=1 𝑟1 (

∑

𝑚=1 𝑟1

∑

𝑚=1 𝑟1

∑

𝑚=1 𝑟1

∑

− 𝑗𝐿𝜋2

2 2

− 𝑗𝐿𝜋2

2 2

) )

𝐿

𝐹01 (𝑚, 𝑗 )𝐹11 (𝑚) ⋅ ∫0 𝑓𝐷1 ′′ (𝑗, 𝜉) sin

𝑖𝜋𝜉 d𝜉, 𝐿

𝐿

𝐹01 (𝑚, 𝑗 )𝑓𝐷 (𝑗)𝐹11 (𝑚) ⋅ ∫0 𝐶𝑜𝑒𝑓 ′′ 1 (𝜉) sin

( ) 𝐿 𝐹01 (𝑚, 𝑗 )𝑓𝐷 (𝑗)Φ𝑚 𝑘𝑚 , ℎ𝑏 ⋅ ∫0 𝐶𝑜𝑒𝑓 ′′ 2 (𝜉) sin 𝐿( 1 𝑗𝜋

𝑖𝜋𝜉 d𝜉, 𝐿

𝑖𝜋𝜉 d𝜉, 𝐿

) ( ) 𝐿 − (−1) 𝐹01 (𝑚, 𝑗 )Φ𝑚 𝑘𝑚 , ℎ𝑏 ⋅ ∫0 𝐶𝑜𝑒𝑓 ′′ 3 (𝜉) sin 𝑗

( ) 𝐿 𝐹01 (𝑚, 𝑗 )Φ𝑚 𝑘𝑚 , ℎ𝑏 ⋅ ∫0 𝑓𝐷1 ′′ (𝑗, 𝜉) sin

𝑚=1 ( ) 𝐿 ∫0 𝑁𝑜𝑛′′ 𝜉; 𝑥0 , 𝑦0 sin 𝑖𝜋𝜉 d𝜉, 𝐿

(B1) 𝑖𝜋𝜉 d𝜉, 𝐿

𝑖𝜋𝜉 d𝜉, 𝐿

𝑛 = 1, 2, ⋯ , 𝑟2 .

References [1] Sakha MS, Shaker HR, Tahavori M. Optimal sensors and actuators placement for large-scale switched systems. Int J Dyn Control 2018:1–10. [2] Fisser M, Badcock RA, Teal PD, Hunze A. Optimizing the sensitivity of palladium based hydrogen sensors. Sens Actuators B Chem 2018;259:10–19. [3] Abbasnejad B, Thorby W, Razmjou A, Jin D, Asadnia M, Warkiani ME. MEMS piezoresistive flow sensors for sleep apnea therapy. Sens Actuators A Phys 2018;279:577–85. [4] Ferrari M, Ferrari V, Guizzetti M, Marioli D, Taroni A. Piezoelectric multifrequency energy converter for power harvesting in autonomous microsystems. Sens Actuators A Phys 2008;142:329–35. [5] Domingo-Roca R, Tiller B, Jackson JC, Windmill JFC. Bio-inspired 3D-printed piezoelectric device for acoustic frequency selection. Sens Actuators A Phys 2018;271:1–8. [6] Buxi D, Redouté JM, Yuce MR. Frequency sensing of medical signals using low-voltage piezoelectric sensors. Sens Actuators A Phys 2014;220:373–81. [7] Park JM, Lee SI, Kwon OY, Choi HS, Lee JH. Comparison of nondestructive microfailure evaluation of fiber-optic Bragg grating and acoustic emission piezoelectric sensors using fragmentation test. Compos Part A Appl Sci Manuf 2003;34:203–16.

369