Nonlinear analysis of axially loaded piezoelectric energy harvesters with flexoelectricity

International Journal of Mechanical Sciences 173 (2020) 105473

Contents lists available at ScienceDirect

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

Nonlinear analysis of axially loaded piezoelectric energy harvesters with flexoelectricity Yunbin Chen a,b, Zhi Yan a,b,∗ a b

Department of Mechanics, Huazhong University of Science and Technology, Wuhan 430074, China Hubei Key Laboratory of Engineering Structural Analysis and Safety Assessment, Luoyu Road 1037, Wuhan 430074, China

a r t i c l e

i n f o

Keywords: Flexoelectricity Axially loaded energy harvester Nonlinear analysis Electromechanical coupling Postbuckling

a b s t r a c t In this work, nonlinear electromechanical models for energy harvesters based on axially preloaded piezoelectric beams incorporating flexoelectric effect are presented. Depending on the amplitude of the applied axial load, the proposed energy harvester could operate in either prebuckling or postbuckling configuration. According to the theory of flexoelectricity and the Hamilton’s principle, the nonlinear electromechanical coupling equations of the proposed energy harvesters under base excitations are derived. For a simply-supported piezoelectric beam, the expression of the static buckling load is analytically determined. Then, for energy harvesters in the prebuckling and postbuckling configurations, we obtain the discrete nonlinear governing equations by employing the Galerkin’s method. These coupling equations are solved numerically by the Runge–Kutta method. Case studies are provided to show the steady-state output voltage and power of the energy harvesters. Results indicate that the frequency response curves show a typical hardening nonlinear behavior in the prebuckling configuration and a softening nonlinear behavior in the postbuckling configuration. Such nonlinear behaviors of energy harvesters imply a wider frequency operation bandwidth of the proposed energy harvesters. We also find that the energy harvesters utilizing flexoelectricity have a better performance than those based on piezoelectricity in the prebuckling configuration. Moreover, we examine the influences of resistive load, mechanical damping coefficient and the amplitude of the base excitation on the performance of the energy harvesters as well as the size effect due to flexoelectricity. It is also interesting to observe that both intrawell and interwell oscillations of the energy harvesters could occur in the postbuckling configuration. This work provides an efficient route to design energy harvesting systems at micro- and nano-scales with enhanced performance.

1. Introduction With the development of advanced manufacturing technology, various electronic devices such as wireless sensors, wearable electronics and micro-electromechanical systems have been designed. Traditionally, these electronic devices require external batteries to provide electrical energy for operation. However, the disadvantages of batteries include low energy density, adverse environmental impacts and the need to either replace or recharge them periodically. As a result, researchers have attempted to harvesting energy from ambient environments to eliminate the need for a battery. Energy harvesting techniques including electrostatic [1], electromagnetic [2] and piezoelectric [3] have been used to convert the ambient mechanical energy to electricity. Among them, piezoelectric energy harvesting technique is most widely adopted. At millimeter scale, Erturk and Inman [4,5] theoretically and experimentally investigated the electromechanical responses of cantilevered unimorph and bimorph piezoelectric energy harvesters. The modeling

∗

was based on the Euler–Bernoulli beam assumption and analytical solutions of the electromechanical responses under the base excitation were obtained and validated by experimental results. At microscale, Xu et al. [6] employed a lumped element model to characterize a PZT thick-film bimorph energy harvester with an integrated silicon proof mass. Kuo et al. [7] fabricated a bimorph MEMS piezoelectric energy harvester with two PZT layers bonded to a stainless steel substrate, experimental results showed that the maximum output power could be hundreds microwatts. Based on ZnO piezoelectric nanowires, a piezoelectric nanogenerator was firstly proposed by Wang and Song [8]. Then, various types of piezoelectric nanogenerators were demonstrated and analyzed [9]. It is evident from these works that piezoelectric energy harvesting technique is popular from macro- to nano-scale and provides a viable route to power electronics requiring minimal amounts of power to function. Piezoelectric energy harvesters proposed in most studies are linear and they operate only in a narrow region of the excitation frequency

Corresponding author at: Department of Mechanics, Huazhong University of Science and Technology, Wuhan 430074, China. E-mail address: [email protected] (Z. Yan).

https://doi.org/10.1016/j.ijmecsci.2020.105473 Received 19 December 2019; Received in revised form 17 January 2020; Accepted 22 January 2020 Available online 24 January 2020 0020-7403/© 2020 Elsevier Ltd. All rights reserved.

Y. Chen and Z. Yan

spectrum, i.e. when the excitation frequency is in the vicinity of the fundamental resonant frequency of the system. In order to obtain a wider frequency operation bandwidth for piezoelectric energy harvesters, researchers have resorted to nonlinear harvesting systems. One typical method is based on cantilevered beams and to use a magnetic mass at the beam’s tip end and extra external magnets to induce nonlinear excitation forces to the beams. For example, Cottone et al. [10] presented a nonlinear piezoelectric oscillator with a magnet attached at the beam tip and an external magnet being placed at an adjustable distance to introduce a repulsive force to the cantilevered beam. Erturk et al. [11] developed a piezomagnetoelastic bimorph generator with two permanent magnets placed near the free end of the cantilevered beam. Tang and Yang [12] experimentally and numerically investigated the performance of a nonlinear piezoelectric energy harvester with a magnetic oscillator. Kim and Seok [13] investigated a multi-stable bimorph cantilever energy harvester using magnetic attraction effect, a bifurcation analysis of the system’s equilibrium was performed and very complex bifurcation scenarios were indicated. Recently, Zhou and Zuo [14] presented a general magnet-coupled piezoelectric tristable energy harvester and analyzed the influence mechanism of asymmetry of potential wells on energy harvesting performance based on the harmonic balance solutions. Jiang et al. [15] employed a 2 degree-of-freedom model to characterize the electromechanical responses of a cantilevered piezoelectric energy harvester with nonlinear magnetic force, experimental work was also provided to verify the model. Another efficient method is to apply an axial load to a beam with fixed or simply-supported ends. For example, Masana and Daqaq [16] investigated the effect of the axial load on the performance of a clampedclamped unimorph piezoelectric energy harvester under transverse base excitations. It was found that a broadband behavior and an enhancement of the power can be achieved due to the compressive axial load. They [17] also compared the performance of a clamped-clamped unimorph piezoelectric energy harvester in the prebuckling (mono-stable) and postbuckling (bi-stable) configurations under harmonic base excitations. Cottone et al. [18,19] experimentally and numerically investigated the performance of a clamped-clamped bimorph, a simply-supported and a clamped-clamped unimorph piezoelectric energy harvesters in the postbuckling configuration under random excitations. Zhu et al. [20] investigated the performance of a simply-supported piezoelectric energy harvester with a magnet-induced compressive axial force. Numerical results revealed that the energy harvester in the postbuckling configuration exhibits a softening nonlinear behavior and can operate in a broader frequency bandwidth. Furthermore, Zhu and Zu [21] developed a magnet-induced buckled-beam piezoelectric generator for broadband vibration-based energy harvesting at low frequencies and small excitations. Numerical simulations indicated strong nonlinearities under harmonic excitations, including snap-through motions, large-amplitude voltage outputs and broad frequency bandwidth. Using the perturbation method, Panyam and Daqaq [22] obtained the approximate solutions of the steady-state periodic electrical responses of both mono-stable and bistable energy harvesters under harmonic base excitations. They found that the bi-stable energy harvesters would produce higher power output than the mono-stable ones for some special cases. More recently, Li et al. [23] proposed a buckled compressive-mode piezoelectric energy harvester under both harmonic and random excitations. The influence of parameters on the performance was discussed via the harmonic balance method and numerical technique. It is also noted in a recent work that lever-bistable energy harvesters utilizing different lever structures could enhance the inter-well dynamic response for improvement of vibration energy harvesting [24]. Tran et al. [25] presented a comprehensive review of various nonlinear techniques to improve the performance of piezoelectric energy harvesters. In recent years, flexoelectricity has attracted significant attention from research communities. The direct flexoelectricity refers to a spontaneous electric polarization generated by a mechanical strain gradient or inhomogeneous strain [26]. When the structural size reduces to micro-

International Journal of Mechanical Sciences 173 (2020) 105473

and nano-scales, flexoelectricity could contribute to the electromechanical coupling of piezoelectric structures, and thus provides a viable route for energy harvesting. Deng et al. [27] proposed a single-layer beam flexoelectric energy harvester under harmonic mechanical excitations. They found that the energy conversion efficiency increases remarkably when the beam thickness reduces from micro- to nano-scale. Wang and Wang [28] presented approximate solutions for the electrical output and optimal load resistance of a unimorph nanoscale piezoelectric energy harvester considering the flexoelectric effect. Yan [29] studied the influences of both flexoelectricity and surface effect on the performance of a single-layer beam nanoscale energy harvester. Recently, Basutkar [30] proposed the analytical solutions for the electrical output and tipvelocity of a bimorph piezoelectric energy harvester with the consideration of flexoelectricity. These flexoelectric energy harvesters provide new insights into the energy harvesting techniques. However, they are linear and thus can operate effectively only in a narrow range of excitation frequencies. By considering geometric nonlinearity, Wang and Wang [31] developed models for both a single-layer and a unimorph cantilevered piezoelectric energy harvester incorporating flexoelectric effect under base excitations. Results showed that the electrical outputs exhibit a nonlinear behavior and a broader bandwidth of operating frequency is achieved. Under force excitations, Dai et al. [32] investigated the electromechanical responses of a single-layer cantilevered beam flexoelectric energy harvester with a mass attached at the beam end. By using a pair of magnets at the tip of a cantilevered plate, Kumar et al. [33] studied the performance of a bi-stable flexoelectric energy harvester using finite element method. In a more recent work, Rojas et al. [34] investigated the influence of flexoelectric and nonlocal elastic effects on the performance of a single-layer beam energy harvester with fixed-hinged boundary condition. In their work, the von-kármán geometric nonlinear strain were considered. Results showed that the power density at macro- to micro-scale is small while is significant at nanoscale due to the size-dependent flexoelectric effect. Besides, Chen et al. [35] studied the performance of a bi-stable plate flexoelectric energy harvester and found that the harvester shows a snap-through behavior, which can lead to the enhancement of the output voltage. It should be mentioned that these broadband flexoelectric energy harvesters are achieved based on geometric nonlinearity of the structure or nonlinear magnetic forces. In this work, we will propose a novel broadband piezoelectric and flexoelectric energy harvester under an applied axial compressive load and base excitations. The energy harvester could operate in either prebuckling configuration or postbuckling configuration, depending on the amplitude of the applied axial load. The modeling of the proposed energy harvester is based on the von-kármán beam assumption and the theory of flexoelectricity. The work is organized as follows. In Section 2, models of piezoelectric energy harvesters under an axial load and base excitations are proposed and the nonlinear electromechanical coupling equations are derived. In Section 3, the solution procedures of static buckling and dynamic electromechanical responses of energy harvesters in prebuckling and postbuckling configurations are presented. Case studies of the performance of energy harvesters are provided and results are discussed in Section 4; and finally, Section 5 concludes this work. 2. Formulations of the problem In this work, we consider a simply-supported piezoelectric beam energy harvester under an axial load and harmonic base excitation. The beam has a length L, width b and a thickness h, respectively. A Cartesian coordinate system (Oxyz) is adopted with x-axis being along the longitudinal direction and z-axis being along the thickness direction of the beam. When the static compressive axial load P is below the first critical 1 , the beam is in a prebuckling configuration, as shown buckling load 𝑃𝑐𝑟 1 , the beam is in a in Fig. 1(a). When the axial load P is larger than 𝑃𝑐𝑟 postbuckling configuration, as shown in Fig. 1(b). The upper and lower

Y. Chen and Z. Yan

International Journal of Mechanical Sciences 173 (2020) 105473

Fig. 1. Schematic of a simply-supported piezoelectric beam energy harvester under an axial load and base excitation; (a) prebuckling configuration, (b) postbuckling configuration.

surfaces of the piezoelectric beam are assumed to be fully covered by conductive electrodes with negligible thickness, which are connected to an external resistive load R. Under the base excitation, a voltage will be generated across the resistive load, which is due to the piezoelectricity and flexoelectricity. Based on the Euler–Bernoulli beam assumption, the displacement fields of the beam can be written as 𝑢(𝑥, 𝑧, 𝑡) = 𝑢0 (𝑥, 𝑡) − 𝑧

𝜕𝑤(𝑥, 𝑡) ; 𝜕𝑥

𝑤(𝑥, 𝑧, 𝑡) = 𝑤(𝑥, 𝑡) + 𝑤𝑏 (𝑡)

(1)

where u0 (x,t) denotes the axial displacement along the central axis of the beam, and w(x, t) represents the transverse displacement relative to the moving base, which has a displacement form of wb (t) = W0 cos (𝜔t) with 𝜔 being the excitation frequency and W0 being the amplitude of the base excitation. Considering the von-kármán geometric nonlinearity [36], the axial strain 𝜀x and the strain gradient 𝜀x,z of the beam can be expressed as ( ) 𝜕𝑢 1 𝜕𝑤 2 𝜕2 𝑤 𝜕2 𝑤 𝜀𝑥 = 0 + −𝑧 ; 𝜀𝑥,𝑧 = − (2) 2 𝜕𝑥 2 𝜕𝑥 𝜕𝑥 𝜕 𝑥2 where the strain gradient 𝜀x,x is neglected here since it is much smaller than 𝜀x,z for a thin beam [37]. The electric field is assumed to exist in the thickness direction only and is expressed as Ez = −V(t)/h. Accordingly, the electric Gibbs free energy density function ub incorporating the flexoelectric effect induced by the strain gradient 𝜀x,z can be expressed as [38] 1 1 1 𝑢𝑏 = − 𝑎33 𝐸𝑧2 + 𝑐11 𝜀2𝑥 − 𝑒31 𝐸𝑧 𝜀𝑥 − 𝑓 𝐸𝑧 𝜀𝑥,𝑧 + 𝑔𝜀2𝑥,𝑧 2 2 2

(3)

where a33 , c11 and e31 are the dielectric, elastic and piezoelectric constants, respectively. f is the direct flexoelectric constant and g denotes the purely non-local elastic effect. Then, the total electric Gibbs free energy Ub in the volume of the beam (Ωb ) is written as 𝑈𝑏 =

{ [ ( ) ] 𝑐 ℎ 𝜕 𝑢0 1 ( 𝜕𝑤 )2 1 𝑉2 − 𝑎33 + 11 + + 𝑒31 𝑉 ∫Ω𝑏 ∫0 2 ℎ 2 𝜕𝑥 2 𝜕𝑥 [ ] ( ) ( 2 )2 } ( ) 3 2 𝜕 𝑢0 1 𝜕𝑤 𝜕 2 𝑤 1 𝑐11 ℎ 𝜕 𝑤 + − 𝑓𝑉 + + 𝑔ℎ d𝑥 (4) 𝜕𝑥 2 𝜕𝑥 2 12 𝜕 𝑥2 𝜕 𝑥2 𝑢 𝑏 d Ω𝑏 = 𝑏

𝐿

Besides, the total kinetic energy of the beam can be written as [( ) ( )2 ] 𝜕 𝑢0 2 1 𝜕𝑤 d𝑤𝑏 𝑇 = 𝜌 + + d Ω𝑏 2 ∫Ω𝑏 𝜕𝑡 𝜕𝑡 d𝑡 [ ] ( ) ( ) 𝐿 ( 𝜕 𝑢 )2 d𝑤𝑏 2 𝜕𝑤 2 𝑚 𝜕𝑤 d𝑤𝑏 0 = + + +2 d𝑥 2 ∫0 𝜕𝑡 𝜕𝑡 d𝑡 𝜕𝑡 d𝑡

(5)

with m = 𝜌bh being the mass per unit length of the beam. In addition, the total work done by the compressive axial load P, damping force and the electric charge output Q(t) can be expressed as ( ) ] 𝐿 [ 𝜕𝑢 𝐿 1 𝜕𝑤 2 𝜕𝑤 0 𝑊 =𝑃 + d𝑥 − 𝑐 𝑤d𝑥 + 𝑄(𝑡)𝑉 (𝑡) (6) ∫0 ∫0 𝑣 𝜕𝑡 𝜕𝑥 2 𝜕𝑥

with cv being the viscous damping coefficient. Substituting Eqs. (4)–(6) into the following Hamilton’s principle, 𝛿

𝑡2

(

∫𝑡1

) 𝑇 − 𝑈𝑏 + 𝑊 d𝑡 = 0

(7)

we obtain

{ [ ( ) ] 𝜕 𝑢0 𝜕 𝛿𝑢0 𝜕𝑤 d𝑤𝑏 𝜕 𝛿𝑤 𝑚 + + ∫𝑡1 ∫0 𝜕𝑡 𝜕𝑡 𝜕𝑡 d𝑡 𝜕𝑡 [ ( ) ]( ) ( ) 2 𝜕 𝑢0 1 𝜕𝑤 𝜕 𝛿𝑢0 𝜕𝑤 𝜕 𝛿𝑤 + 𝑃 − 𝐸 𝐴𝑒𝑓 𝑓 + − 𝑒31 𝑏𝑉 + 𝜕𝑥 2 𝜕𝑥 𝜕𝑥 𝜕𝑥 𝜕𝑥 ( ) 2 𝜕2 𝑤 𝜕 𝛿𝑤 𝜕𝑤 − 𝐸 𝐼𝑒𝑓 𝑓 − 𝑓 𝑏𝑉 − 𝑐𝑣 𝛿𝑤 𝜕𝑡 𝜕 𝑥2 𝜕 𝑥2 [ ( ) ( )] } 𝜕 𝑢0 1 ( 𝜕𝑤 )2 𝜕2 𝑤 𝑏𝑉 − 𝑒31 𝑏 + − 𝑓𝑏 + 𝑎33 𝛿𝑉 (𝑡) d𝑥 𝜕𝑥 2 𝜕𝑥 ℎ 𝜕 𝑥2 𝑡2

d𝑡

𝑡2

+

∫𝑡1

𝐿

𝑄(𝑡)𝛿𝑉 (𝑡)d𝑡 = 0

(8)

with EAeff = c11 bh and EIeff = c11 bh3 /12 + gbh being the effective tensile and bending rigidity of the piezoelectric beam, respectively. Integrating by parts and using the relation i(t) = V(t)/R = dQ(t)/dt, the nonlinear electromechanical coupling equations can be obtained from Eq. (8) as [ ( ) 𝜕 2 𝑢0 𝜕 𝑢0 1 ( 𝜕𝑤 )2 𝜕 𝑚 + 𝑃 − 𝐸 𝐴𝑒𝑓 𝑓 + 𝜕𝑥 𝜕𝑥 2 𝜕𝑥 𝜕 𝑡2 ] −𝑒31 𝑏𝑉 [𝐻(𝑥) − 𝐻(𝑥 − 𝐿)] = 0 (9) ( 𝑚

) 2 𝜕 2 𝑤 d 𝑤𝑏 𝜕𝑤 𝜕 4 𝑤 𝜕 2 [𝐻(𝑥) − 𝐻(𝑥 − 𝐿)] + + 𝑐𝑣 + 𝐸 𝐼𝑒𝑓 𝑓 − 𝑓 𝑏𝑉 𝜕𝑡 𝜕 𝑡2 d𝑡2 𝜕 𝑥4 𝜕 𝑥2 [( ( ) ) ] ( ) 𝜕 𝑢0 1 𝜕𝑤 2 𝜕 𝜕𝑤 + 𝑃 − 𝐸 𝐴𝑒𝑓 𝑓 + − 𝑒31 𝑏𝑉 =0 (10) 𝜕𝑥 𝜕𝑥 2 𝜕𝑥 𝜕𝑥

𝐶𝑝

( ( ) )] 𝐿[ d𝑉 (𝑡) 𝑉 (𝑡) 𝜕3 𝑤 𝜕 𝜕 𝑢0 1 𝜕𝑤 2 + + 𝑓𝑏 − 𝑒31 𝑏 + d𝑥 = 0 ∫0 d𝑡 𝑅 𝜕𝑡 𝜕𝑥 2 𝜕𝑥 𝜕 𝑥2 𝜕𝑡 (11)

with Cp = a33 bL/h being the capacitance of the piezoelectric energy harvester. H(x) is the Heaviside step function and H(x) − H(x − L) is multiplied to make sure that the electromechanical coupling term can survive after a spatial derivative [39]. Besides, the corresponding boundary conditions are { ( ) } 𝜕 𝑢0 1 ( 𝜕𝑤 )2 𝑃 − 𝐸 𝐴𝑒𝑓 𝑓 + − 𝑒31 𝑏𝑉 [𝐻 (𝑥) − 𝐻 (𝑥 − 𝐿)] 𝛿𝑢0 = 0 𝜕𝑥 2 𝜕𝑥 (12) {

[ ] 𝜕 𝜕2 𝑤 𝐸 𝐼𝑒𝑓 𝑓 − [𝐻(𝑥) − 𝐻(𝑥 − 𝐿)]𝑓 𝑏𝑉 𝜕𝑥 𝜕 𝑥2 ( ( ) ) } 𝜕 𝑢0 1 ( 𝜕𝑤 )2 𝜕𝑤 + 𝑃 − 𝐸 𝐴𝑒𝑓 𝑓 + − 𝑒31 𝑏𝑉 𝛿𝑤 = 0 𝜕𝑥 2 𝜕𝑥 𝜕𝑥

(13)

Y. Chen and Z. Yan

[

International Journal of Mechanical Sciences 173 (2020) 105473

] 𝜕2 𝑤 𝜕 𝛿𝑤 𝐸 𝐼𝑒𝑓 𝑓 − [𝐻(𝑥) − 𝐻(𝑥 − 𝐿)]𝑓 𝑏𝑉 =0 𝜕𝑥 𝜕 𝑥2

(14)

at x = 0 and x = L. For a slender beam, the longitudinal inertia term can be ignored, 𝜕2 𝑢 𝑚 𝜕 𝑡20

i.e. ≈ 0. As a result, we can obtain the following expression from Eq. (9) as 𝜕 𝑢0 1 ( 𝜕𝑤 )2 + = 𝑐1 (𝑡) (15) 𝜕𝑥 2 𝜕𝑥 Using the boundary condition in Eq. (12), c1 (t) is determined as 1 2 𝐿 ∫0

𝑐1 =

𝐿

(

𝜕𝑤 𝜕𝑥

)2

d𝑥

(16)

By introducing the following nondimensional parameters √ √ 𝐸 𝐼𝑒𝑓 𝑓 𝑥 𝑤 𝑚𝐿4 𝑋= , 𝑊 = , 𝑇 =𝑡 , Ω=𝜔 𝐿 𝑟 𝐸 𝐼𝑒𝑓 𝑓 𝑚𝐿4

1

∫0

[ Θ 𝑊̇ ′′ − Θ𝑝 𝑓

1

(17)

∫0

] 𝑊 𝑊̇ ′ d𝑋 d𝑋 = 0 ′

(19)

with 𝑐𝑣 𝐿2

⌢ 𝑒 𝑏𝐿2 𝑓 𝑓 𝑏𝐿2 𝑃 𝐿2 , 𝜗𝑝 = 31 , 𝜗 = , )1∕2 , 𝑃 = 𝐸 𝐼 𝐸 𝐼𝑒𝑓 𝑓 𝑟𝐸 𝐼𝑒𝑓 𝑓 𝑒𝑓 𝑓 𝑚𝐸 𝐼𝑒𝑓 𝑓 √ 𝑚𝜔2 𝑊0 𝐿4 𝑒 𝑏𝑟2 𝑓 𝑏𝑟 𝑚𝐿4 1 𝐹 = , 𝛼= , Θ𝑓 = , Θ𝑝 = 31 𝑟𝐸 𝐼𝑒𝑓 𝑓 𝐸 𝐼𝑒𝑓 𝑓 𝑅𝐶𝑝 𝐿𝐶𝑝 𝐿𝐶𝑝

𝑐 = (

(20)

⌢

the following analysis, we use the notation P instead of 𝑃 for convenience. For a simply-supported beam considered in the current work, the associated nondimensional boundary conditions can be obtained from Eqs. (13) and (14) as (21)

under the short circuit condition (i.e. R → 0, V → 0).

where 𝜆2 = P − Γ. Solving the above equation and using the boundary conditions in Eq. (21), the static buckling configuration of the simplysupported beam can be expressed as 𝑊𝑠 (𝑋) = 𝐶𝑠 sin(𝜆𝑋)

where 𝜆 = ns 𝜋 with ns = 1, 2...N being an integer representing the number of static buckling mode. Using the relation 𝜆2 = 𝑃 − Γ = 𝑃 −

1

1 2 𝑊𝑠 ′ d𝑋 2 ∫0

(26)

√ 𝐶𝑠 = ±2

𝑃 −1 𝜆2

(27)

which requires the axial load greater than the critical buckling load, i.e. 𝑃 ≥ 𝑃𝑐𝑟 = 𝑛2𝑠 𝜋 2 . Thus, the static buckling configurations corresponding to a given axial load P under simply-supported boundary condition is determined. 3.2. Electromechanical responses of the energy harvester in the prebuckling configuration When the applied axial load P is less than the first critical buckling 1 = 𝜋 2 , the beam is in the prebuckling configuration. To load, i.e. 𝑃 < 𝑃𝑐𝑟 investigate the electromechanical coupling behaviors of the energy harvester, the natural frequency and the trial function of the piezoelectric beam are firstly determined. We assume that (28)

where Ω is the undamped natural frequency and 𝜙(X) is the corresponding mode shape. Substituting Eq. (28) into (18) and ignoring the damping, external force and nonlinear terms, we obtain 𝜙(𝑖𝑣) + 𝑃 𝜙′′ − Ω2 𝜙 = 0

(29)

the general solution is given by 𝜙(𝑋) = 𝐶1 cos(𝑠1 𝑋) + 𝐶2 sin(𝑠1 𝑋) + 𝐶3 𝑒𝑠2 𝑋 + 𝐶4 𝑒−𝑠2 𝑋

(30)

with

√ √ √( √( )1 )1 √ √ 2 2 √ 𝑃2 √ 𝑃2 𝑃 𝑃 𝑠1 = + Ω2 + ; 𝑠2 = + Ω2 − 4 2 4 2

(31)

Using the boundary conditions in Eq. (21) yields the following algebraic equations

𝐶1 cos(𝑠1 ) + 𝐶2 sin(𝑠1 ) + 𝐶3 𝑒𝑠2 + 𝐶4 𝑒−𝑠2 = 0 −𝐶1 𝑠21 + 𝐶3 𝑠22 + 𝐶4 𝑠22 = 0

3.1. Static buckling analysis

−𝐶1 𝑠21 cos(𝑠1 ) − 𝐶2 𝑠21 sin(𝑠1 ) + 𝐶3 𝑠22 𝑒𝑠2 + 𝐶4 𝑠22 𝑒−𝑠2 = 0

For the short circuit condition (i.e. R → 0, V → 0), the static buckling problem can be obtained by dropping the time-related, damping and external forcing terms [40], then Eq. (18) is simplified as ( ) 1 1 2 𝑊𝑠 (𝑖𝑣) + 𝑃 − 𝑊𝑠 ′ d𝑋 𝑊𝑠 ′′ = 0 (22) 2 ∫0 where Ws denotes the static buckled configuration. Since the integral term in the above equation is a constant, we let

sin(𝑠1 ) = 0

(33)

which means s1 = n𝜋 with n being an integer corresponding to the vibrational mode shape. Substituting the value of s1 = n𝜋 into Eq. (32) and using the following normalized condition of the vibrational mode shapes [41] 1

(23)

(32)

Accordingly, the determinant of the coefficient matrix should be zero and we can obtain

1

1 2 𝑊𝑠 ′ d𝑋 2 ∫0

(25)

𝐶1 + 𝐶3 + 𝐶4 = 0

3. Solution procedure

Γ=

(24)

𝑊 (𝑋, 𝑇 ) = 𝜙(𝑋)e𝑖Ω𝑇

The prime in Eqs. (18) and (19) stands for a derivation with respect to spatial coordinate X and the overdot indicates the derivative with respect to time T. Coefficients ϑp and Θp denote the electromechanical coupling terms due to piezoelectricity, while ϑf and Θf denote the electromechanical coupling terms due to flexoelectricity. In

𝑊 = 𝑊 ′′ = 0 at 𝑋 =0 and 𝑋 =1

𝑊𝑠 (𝑖𝑣) + 𝜆2 𝑊𝑠 ′′ = 0

we can determine the amplitude of the static buckling configuration as

√ where 𝑟 = 𝐸 𝐼𝑒𝑓 𝑓 ∕𝐸 𝐴𝑒𝑓 𝑓 . The nondimensional forms of the electromechanical coupling equations with the consideration of Eq. (15) are then expressed as ( ) 1 ⌢ 1 (𝑖𝑣) ′2 ̈ ̇ 𝑊 + 𝑐𝑊 + 𝑊 + 𝑃− 𝑊 d𝑋 𝑊 ′′ 2 ∫0 [ [ ]] − 𝜗𝑝 𝑊 ′′ + 𝜗𝑓 𝐻 ′′ (𝑋) − 𝐻 ′′ (𝑋 − 1) 𝑉 (𝑡) = 𝐹 cos(Ω𝑇 ) (18)

𝑉̇ + 𝛼𝑉 +

Then, Eq. (22) reduces to a fourth-order ordinary differential equation with constant coefficients

∫0

𝜙𝑖 𝜙𝑗 d𝑋 = 𝛿𝑖𝑗

(34)

Y. Chen and Z. Yan

International Journal of Mechanical Sciences 173 (2020) 105473

where 𝛿 ij is the Kronecker delta, the nth normalized mode shape can be obtained as √ 𝜙𝑛 (𝑋) = 2 sin(𝑛𝜋𝑋) (35) Besides, using the relation defined in Eq. (31), the corresponding nth natural frequency is obtained as √ Ω𝑛 = 𝑛𝜋 𝑛2 𝜋 2 − 𝑃 (36) which requires P ≤ n2 𝜋 2 . Adopting the Galerkin’s method, the nondimensional transverse displacement of the unbuckled beam can be expressed as 𝑊 (𝑋, 𝑇 ) =

𝑁 ∑ 𝑗=1

𝑞𝑗 (𝑇 )𝜙𝑗 (𝑋)

(37)

where N is the number of modes used in the series discretization and qj (T) is the unknown generalized coordinate. Substituting Eq. (37) into the electromechanical coupling Eqs. (18) and (19), multiplying by 𝜙j and integrating the equations from X = 0 to 1, the following discrete Euler– Lagrange equations of energy harvester in the prebuckling configuration are obtained as ( ) 𝑞̈𝑖 + 𝜂𝑖𝑗 𝑞̇ 𝑗 + 𝜅𝑖𝑗 − 𝜁𝑖𝑗𝑝 𝑉 𝑞𝑗 − 𝐵𝑖𝑗𝑘𝑙 𝑞𝑗 𝑞𝑘 𝑞𝑙 − 𝜍𝑖𝑓 𝑉 = 𝑓𝑖 cos(Ω𝑇 ) (38) 𝑉̇ + 𝛼𝑉 +

𝛽𝑖𝑓 𝑞̇ 𝑖

−

𝛾𝑖𝑗𝑝 𝑞𝑖 𝑞̇ 𝑗

=0

(39)

with 1

𝜂𝑖𝑗 = 𝑐 𝜅𝑖𝑗 =

∫0 1

∫0

𝑓𝑖 = 𝛽𝑖𝑓

[ ] 𝜙𝑖 𝜙(𝑗𝑖𝑣) + 𝑃 𝜙𝑗 ′′ d𝑋 = Ω2𝑗 𝛿𝑖𝑗 ( ) 1

1

𝜁𝑖𝑗𝑝 = 𝜗𝑝 =𝜗

(40a)

∫0 1

𝑓

∫0

1

𝜙𝑖 𝜙𝑗 ′′ d𝑋 [

(40d)

𝜙𝑖 𝐹 d𝑋

∫0

[ ] = Θ 𝜙𝑖 ′ (1) − 𝜙𝑖 ′ (0) 𝑓

𝛾𝑖𝑗𝑝 = Θ𝑝

1

∫0

𝜙𝑖 ′ 𝜙𝑗 ′ d𝑋

𝑊𝑠 ′ 𝑊𝑑 ′ d𝑋

(44)

Similarly, we let 𝑊𝑑 (𝑋, 𝑇 ) = 𝜙(𝑋)e𝑖Ω𝑇

(45)

where Ω is the undamped natural frequency and 𝜙(X) is the corresponding mode shape. Substituting Eq. (45) into Eq. (44), the following nonhomogeneous fourth-order ordinary differential equation is obtained as 𝜙(𝑖𝑣) + 𝜆2 𝜙′′ − Ω2 𝜙 = 𝑊𝑠 ′′

1

∫0

𝑊𝑠 ′ 𝜙′ d𝑋

(46)

The solution of Eq. (46) can be expressed as (47)

in which, the homogeneous solution is ( ) ( ) ( ) ( ) 𝜙ℎ (𝑋 ) = 𝑑1 sin 𝑆1 𝑋 + 𝑑2 cos 𝑆1 𝑋 + 𝑑3 sinh 𝑆2 𝑋 + 𝑑4 cosh 𝑆2 𝑋 (48)

𝑆1 =

[

𝜆2 +

√

𝜆4 + 4Ω2 2

]1∕2

[ ; 𝑆2 =

−𝜆2 +

]1∕2 √ 𝜆4 + 4Ω2 2

and a particular solution has the form of [40]

(40g)

𝜙𝑝 (𝑋) = 𝑑5 𝑊𝑠 ′′

1

(41)

3.3. Electromechanical responses of the energy harvester in the postbuckling configuration When the applied axial load P is greater than the first critical buckling load, the piezoelectric beam will buckle. To investigate the electromechanical responses of the energy harvester in the postbuckling configuration, we let (42)

(49)

(50)

Substituting Eqs. (47) into Eq. (46) and using Eq. (24) yields ( )

(40h)

where G(q) represents the geometric nonlinear term and F(T) denotes the external force due to the base excitation. To this end, the generalized coordinate qj (T) and the output voltage V of the energy harvester in the prebuckling configuration can be numerically solved by the Runge– Kutta method, which is implemented by using Matlab function ode45 in the present work.

𝑊 (𝑋, 𝑇 ) = 𝑊𝑠 (𝑋) + 𝑊𝑑 (𝑋, 𝑇 )

1

∫0

(40f)

By introducing a new state vector 𝐲 = {𝑞, 𝑞̇ , 𝑉 }, the discrete Euler– Lagrange equations can be rewritten into the following implicit firstorder form 𝐲̇ = 𝐀𝐲 + 𝐆(𝑞) + 𝐅(𝑇 )

𝑊̈ 𝑑 + 𝑊𝑑(𝑖𝑣) + 𝜆2 𝑊𝑑 ′′ = 𝑊𝑠 ′′

with (40e)

(43)

In order to obtain the trial functions, we investigate the linear free vibration of the beam in the buckled configuration first. The governing equation is simplified by dropping the damping, external force and nonlinear terms in Eq. (43), and we obtain the following equation under the short circuit condition (i.e. R → 0, V → 0)

𝜙(𝑋) = 𝜙ℎ (𝑋) + 𝜙𝑝 (𝑋) (40c)

]

[ ] d𝛿(𝑋) d𝛿(𝑋 − 1) 𝜙𝑖 − d𝑋 = 𝜗𝑓 𝜙𝑖 ′ (1) − 𝜙𝑖 ′ (0) d𝑋 d𝑋

𝑊̈ 𝑑 + 𝑐 𝑊̇ 𝑑 + 𝑊𝑑(𝑖𝑣) + 𝜆2 𝑊𝑑 ′′ [( ) 1 1 1 ′ ′ ′2 − 𝑊𝑠 𝑊𝑑 d𝑋 + 𝑊𝑑 d𝑋 𝑊𝑠 ′′ ∫0 2 ∫0 ( ) ] 1 1 1 ′ ′ ′2 ′′ + 𝑊𝑠 𝑊𝑑 d𝑋 + 𝑊𝑑 d𝑋 𝑊𝑑 ∫0 2 ∫0 [ 𝑝 ( ′′ ) [ ]] − 𝜗 𝑊𝑠 + 𝑊𝑑 ′′ + 𝜗𝑓 𝐻 ′′ (𝑋) − 𝐻 ′′ (𝑋 − 1) 𝑉 = 𝐹 cos(Ω𝑇 )

(40b)

1

1 𝜙 𝜙 ′ 𝜙 ′ d𝑋 𝜙𝑙 ′′ d𝑋 2 ∫0 𝑖 ∫0 𝑗 𝑘

𝐵𝑖𝑗𝑘𝑙 =

𝜍𝑖𝑓

𝜙𝑖 𝜙𝑗 d𝑋 = 𝑐 𝛿𝑖𝑗

where Ws (X) is the static buckling configuration and Wd (X,T) is the dynamic disturbance around the static buckling configuration. Substituting the above equation into Eq. (18) and using Eq. (22), we obtain

∫0

𝑊𝑠 ′ 𝜙ℎ ′ d𝑋 + 𝑑5 Ω2 +

1

∫0

𝑊𝑠 ′ 𝑊𝑠 ′′′ d𝑋

=0

(51)

Using the boundary conditions in Eq. (21), the following algebraic equations are obtained as 𝑑2 + 𝑑4 = 0 −𝑑2 𝑆12 + 𝑑4 𝑆22 = 0 ( ) ( ) ( ) ( ) 𝑑1 sin 𝑆1 + 𝑑2 cos 𝑆1 + 𝑑3 sinh 𝑆2 + 𝑑4 cosh 𝑆2 = 0 ( ) ( ) ( ) ( ) −𝑑1 𝑆12 sin 𝑆1 − 𝑑2 𝑆12 cos 𝑆1 + 𝑑3 𝑆22 sinh 𝑆2 + 𝑑4 𝑆22 cosh 𝑆2 = 0 (52) where d1 -d5 can be determined from the algebraic Eqs. (51) and (52) with the consideration of Eqs. (25) and (48). To obtain nontrivial solutions, the following characteristic equation should be satisfied ( ) 𝐶 2 𝜆4 ( 2 )2 ( ) ( ) Ω2 − 𝑠 (53) 𝑆1 + 𝑆22 sin 𝑆1 sinh 𝑆2 = 0 2 If Ω2 −

𝐶𝑠2 𝜆4 2

=0

(54)

Y. Chen and Z. Yan

International Journal of Mechanical Sciences 173 (2020) 105473

the natural frequency Ω is determined and the nth vibrational mode shape with the consideration of the normalized condition of Eq. (34) can be written as √ ( ) 𝜙𝑛 (𝑋 ) = 2sin 𝑛𝑠 𝜋𝑋 , 𝑛 = 𝑛𝑠 (55) Else if Ω2 −

𝐶𝑠2 𝜆4 2

≠ 0, we have

sin(𝑆1 ) = 0 i.e. 𝑆1 = 𝑛𝜋

(56)

or ( ) sinh 𝑆2 = 0 i.e. 𝑆2 = in𝜋

(57)

Substituting Eqs. (56) and (57) into Eq. (49), we obtain the nth nondimensional natural frequency as ( ) Ω2𝑛 = 𝑛2 𝜋 4 𝑛2 − 𝑛2𝑠 , 𝑛 ≠ 𝑛𝑠 (58) when n > ns , Ω2𝑛 is positive and Eq. (56) holds; when n < ns , Ω2𝑛 is negative and Eq. (57) holds. Accordingly, the normalized mode shape can be derived as √ 𝜙𝑛 (𝑋) = 2 sin(𝑛𝜋𝑋), 𝑛 ≠ 𝑛𝑠 (59)

Fig. 2. Variation of the dimensionless fundamental natural frequency with the axial load in the prebuckling and postbuckling domains.

Then, using the Galerkin’s method, the nondimensional transverse displacement of the buckled beam can be expressed as

4. Results and discussions

𝑊 (𝑋, 𝑇 ) = 𝑊𝑠 (𝑋) + 𝑊𝑑 (𝑋, 𝑇 ) = 𝑊𝑠 (𝑋) +

𝑁 ∑ 𝑗=1

𝑞𝑗 (𝑇 )𝜙𝑗 (𝑋)

(60)

Substituting Eq. (60) into the electromechanical coupling Eqs. (18) and (19), multiplying by 𝜙j and integrating the equations from X = 0 to 1, the following discrete Euler–Lagrange equations of the energy harvester in the postbuckling configuration are obtained as ( ) ⌢ 𝑞̈𝑖 + 𝜂𝑖𝑗 𝑞̇ 𝑗 + 𝜅 𝑖𝑗 𝑞𝑗 − 𝐴𝑖𝑗𝑘 𝑞𝑗 𝑞𝑘 − 𝐵𝑖𝑗𝑘𝑙 𝑞𝑗 𝑞𝑘 𝑞𝑙 − 𝜁𝑖𝑗𝑝 𝑞𝑗 + 𝜍𝑖𝑝 + 𝜍𝑖𝑓 𝑉 = 𝑓𝑖 cos(Ω𝑇 )

(61)

( ) 𝑉̇ + 𝛼𝑉 + 𝛽𝑖𝑓 − 𝛽𝑖𝑝 𝑞̇ 𝑖 − 𝛾𝑖𝑗𝑝 𝑞𝑖 𝑞̇ 𝑗 = 0 with 𝜅 𝑖𝑗 =

[

1

⌢

∫0

𝐴𝑖𝑗𝑘 =

(

𝜙𝑖 𝜙(𝑗𝑖𝑣) + 𝜆2 𝜙𝑗 ′′ −

1

∫0

(62) )

]

𝑊𝑠 ′ 𝜙𝑗 ′ d𝑋 𝑊𝑠 ′′ d𝑋 = Ω2𝑗 𝛿𝑖𝑗

(63a)

[( ) ] 1 1 1 𝜙𝑖 𝜙𝑗 ′ 𝜙𝑘 ′ d𝑋 𝑊𝑠 ′′ d𝑋 ∫0 2 ∫0 [( ) ] 1

+

∫0 1

𝜍𝑖𝑝 = 𝜗𝑝

∫0

𝛽𝑖𝑝 = Θ𝑝

∫0

1

𝜙𝑖

1

∫0

𝑊𝑠 ′ 𝜙𝑗 ′ d𝑋 𝜙𝑘 ′′ d𝑋

(63b)

𝜙𝑖 𝑊𝑠 ′′ d𝑋

(63c)

𝑊𝑠 ′ 𝜙𝑖 ′ d𝑋

(63d)

The other coefficients are the same as those in Eqs. (40a) and (40c)– (40h). It is seen that an additional quadratic nonlinear term Aijk qj qk is included in comparison to that in Eq. (38). Besides, two electromechanical coupling terms 𝜍𝑖𝑝 𝑉 and 𝛽𝑖𝑝 𝑞̇ 𝑖 due to the piezoelectric effect are included, which are related to the static buckling configuration Ws . Therefore, it is expected that the electromechanical responses in the postbuckling configuration could be very different from those in the prebuckling configuration. By using the Runge–Kutta method, the problem can be numerically solved.

In this section, case studies are provided to show the performance of the proposed axially loaded energy harvesters with the consideration of flexoelectric effect. The geometric dimension of the piezoelectric beam is fixed to be L:b:h = 100:10:1, and the structural material is chosen to be lead zirconate titanate (i.e. PZT-5H). The material properties of PZT-5H are c11 = 126 GPa, e31 = −6.5 C/m2 and a33 = 1.3 × 10−8 C/(V • m). The mass density is 𝜌 = 7500 kg/m3 . Besides, the direct flexoelectric constant is taken as f = −1.0 × 10−7 C/m [42] and the non-local elastic constant is assumed to be 𝑔 = 𝑐11 𝑙02 with the material length scale l0 = 0.1h [32]. In addition, the nondimensional damping coefficient c in Eq. (18) is set to be c = 0.05 without specifically indicated. Firstly, the variation of the dimensionless fundamental natural frequency with the compressive axial load in both prebuckling and postbuckling domains is shown in Fig. 2. The short circuit condition is adopted. In the prebuckling domain, it is seen that the natural frequency decreases with the increase of the axial load until it reaches zero at P = 𝜋 2 , which corresponds to the first critical buckling load. In the postbuckling domain, the natural frequency increases with the increase of the axial load. These phenomenon are consistent with those observed in the experimental work [17]. Nextly, we will analyze the performance of the energy harvester in the prebuckling configuration, and the performance of the energy harvester in the postbuckling configuration around the first static buckling mode. Since the amplitude of vibration is small and the first mode vibration is dominant, N = 1 is taken in the Galerkin’s discretization, this one-mode approximation has also been adopted in references [17,18]. The variations of the steady-state solution of the maximum voltage with varying nondimensional excitation frequency for different load resistances are plotted in Fig. 3. We choose the corresponding axial loads to make the nondimensional fundamental natural frequencies of the piezoelectric beam in the prebuckling and postbuckling configurations being Ω = 2, as shown in Fig. 2. The amplitude of the base excitation is set to be W0 /h = 0.25%, and the thickness of the beam is chosen to be h = 300 nm. In the current study, both piezoelectric and flexoelectric effects will contribute to the electromechanical coupling. Three groups of results are presented in Fig. 3 to examine the electromechanical couplings contributed by the piezoelectric effect and flexoelectric effect, respectively. Fig. 3(a) and (b) consider the piezoelectric effect (PE) only and we set the flexoelectric coefficient f = 0. Fig. 3(c) and (d) consider the flexoelectric effect (FE) only and we set the piezoelectric coefficient e31 = 0. Fig. 3(e) and (f) consider both the piezoelectric and flexoelectric

Y. Chen and Z. Yan

International Journal of Mechanical Sciences 173 (2020) 105473

Fig. 3. Frequency response curves for the maximum voltage of the energy harvesters with different load resistances in the prebuckling configuration (left column) and postbuckling configuration (right column): (a) and (b) with piezoelectric effect (PE) only; (c) and (d) with flexoelectric effect (FE) only; (e) and (f) with both PE and FE.

effects (PE and FE). It is observed that the maximum voltage output increases with the increase of the load resistance (e.g. from R = 0.2 MΩ to R = 2 MΩ) for all the cases studied, and the voltage will reach to a converged value when R is large enough. We can also see from Fig. 3(a) and (b) that the maximum voltage of the energy harvester in the postbuckling configuration due to the piezoelectric effect (about 0.281 mV) is larger than that in the prebuckling configuration (about 0.237 mV). This phenomenon is due to an additional quadratic nonlinear term Aijk qj qk and two electromechanical terms 𝜍𝑖𝑝 𝑉 and 𝛽𝑖𝑝 𝑞̇ 𝑖 as seen in Eqs. (61) and (62) in the postbuckling configuration as compared with that in the pre-

buckling configuration. Through numerical simulations, it is found that the inclusion of Aijk qj qk will decrease the output voltage, while the additional electromechanical terms will increase the voltage. As a result, the voltage response is slightly enhanced in the postbuckling configuration. In contrast, the maximum voltage of the energy harvester due to the flexoelectric effect in the prebuckling configuration (about 0.527 mV) is more than two times that of the value in the postbuckling configuration (about 0.233 mV), as indicated in Fig. 3(c) and (d), which is due to the incorporation of the term Aijk qj qk . When both piezoelectric and flexoelectric effects are considered, voltage output of energy harvesters

Y. Chen and Z. Yan

International Journal of Mechanical Sciences 173 (2020) 105473

Fig. 4. Variations of the power output of the energy harvesters in the (a) prebuckling and (b) postbuckling configurations with the load resistance.

Fig. 5. Frequency response curves for the maximum voltage of the energy harvesters in the (a) prebuckling and (b) postbuckling configurations with different damping coefficients.

Fig. 6. Frequency response curves for the maximum voltage of the energy harvesters in the (a) prebuckling and (b) postbuckling configurations with different amplitudes of the base excitation. The damping coefficient is set as c = 0.05.

in either case is enhanced accordingly. As compared to that shown in Fig. 3(a), a wider frequency operation bandwidth is achieved in the later cases. In addition, the frequency response curves of the energy harvester in the prebuckling configuration display a typical hardening nonlinear behavior, i.e. the resonant frequency shifts to the right, as seen from

Fig. 3(a), (c) and (e). While the frequency response curves of the energy harvester in the postbuckling configuration show a softening nonlinear behavior, i.e. the resonant frequency shifts to the left, as seen from Fig. 3(b), (d) and (f). This switch from a hardening nonlinear behavior to a softening one is due to the change of the nonlinearity of the harvesting system, i.e. an additional quadratic nonlinear term Aijk qj qk is included in

Y. Chen and Z. Yan

International Journal of Mechanical Sciences 173 (2020) 105473

Fig. 7. Time histories of the mid-point displacement of the energy harvester in the postbuckling configuration (left column) and corresponding phase diagrams (right column) under different amplitudes of the base excitations: (a) and (b) Intrawell oscillation when W0 /h = 0.25%; (c) and (d) Chaotic oscillation when W0 /h = 1.0%; (e) and (f) Interwell oscillation when W0 /h = 2.5%.

the postbuckling configuration in comparison to the prebuckling configuration. Such a hardening nonlinear behavior in the prebuckling configuration and a softening behavior in the postbuckling configuration have also been observed in the reference [17]. Fig. 4(a) and (b) show the variations of the power output of the energy harvesters in the prebuckling and postbuckling configurations with the load resistance, respectively. The power output is defined as 2 ∕𝑅, where V 𝑃 = 𝑉max max is the peak value of the voltage output for a fixed load resistance R during both forward and backward path sweeps

of the nondimensional excitation frequency between Ω = 1 and Ω = 3. Usually, Vmax is obtained within the resonant frequency region around Ω = 2. It is observed that the power output increases significantly with the increase of the load resistance at first, and decreases with the further increase of the load resistance for all the cases studied. Due to the competition mechanism of the voltage and the current (i = V/R, which decreases with the increase of R), there exists an optimal load resistance corresponding to the peak power output. It should be pointed out that the optimal load resistance Ropt is obtained within the reso-

Y. Chen and Z. Yan

nant frequency region. When flexoelectric effect is considered for energy harvesters in both configurations, the optimal load resistance increases slightly, e.g. from 0.15 MΩ to 0.25 MΩ in the prebuckling configuration and from 0.40 MΩ to 0.45 MΩ in the postbuckling configuration. It is also seen that comparing to the piezoelectric effect, the maximum power output due to the flexoelectric effect is much larger in the prebuckling configuration, while is smaller in the postbuckling configuration. The maximum power output is significantly enhanced when both effects are considered. When both piezoelectric and flexoelectric effects are considered, the influences of the mechanical damping coefficient on the maximum output voltage are investigated, as shown in Fig. 5(a) and (b), respectively. The load resistance is set to be R = 2 MΩ. It is observed that with the decrease of the damping coefficient, the maximum voltage increases and the resonant frequency of the energy harvester in the prebuckling configuration shifts to right while shifts to left in the postbuckling configuration. The influence of the damping coefficient on the voltage output is significant only when the excitation frequency is within the frequency resonant region. It is also noted that the frequency response curve of the energy harvester in the prebuckling configuration is linear while is nonlinear in the postbuckling configuration when the damping coefficient is relatively large (e.g. c = 0.10), which indicates that a smaller damping coefficient (e.g. c = 0.05) is required to trigger the nonlinear behavior for the energy harvester in the prebuckling configuration. Fig. 6(a) and (b) show the influences of the amplitude of base excitation on the maximum voltage output of the proposed harvesters. We can see that the maximum voltages are enhanced with the increase of the amplitude of the base excitation in a wide range of the nondimensional excitation frequencies, i.e. 1.7 ≤ Ω ≤ 2.7 in the prebuckling configuration and 1.5 ≤ Ω ≤ 3 in the postbuckling configuration, not just within the resonant frequency region as shown in Fig. 5. It is also seen that when the amplitude is relatively small (e.g. W0 /h = 0.15%), the frequency response curve of the energy harvester in the postbuckling configuration shows a softening nonlinear behavior while is linear in the prebuckling configuration, which means that in the postbuckling configuration, a nonlinear electromechanical response with a wider frequency operation range can be triggered more easily. It should be mentioned that the above dynamic response of the energy harvester in the postbuckling configuration corresponds to the intrawell oscillation. When the amplitude of base excitation increases, interwell oscillation may occur. Fig. 7 plots the steady-state time histories and corresponding phase diagrams of the energy harvester under different amplitudes of the base excitations, the nondimensional excitation frequency is chosen to be Ω = 1.6. As seen from the figure, when the amplitude of the base excitation is small (i.e. W0 /h = 0.25%), the energy harvester periodically vibrates around the static buckling configuration and the phase diagram is a simple circle for the case. This vibrational mode is known as the intrawell oscillation. As the amplitude of base excitation increases to W0 /h = 1.0%, the harvester experiences aperiodic or chaotic vibrations, as seen from Fig. 7(c) and (d). It is further observed from Fig. 7(e) and (f) that periodic interwell oscillation occurs at W0 /h = 2.5% and the energy harvester experiences a snap-through from one stable equilibrium to the other stable equilibrium. These three vibrational modes belong to typical dynamic behaviors of a bi-stable oscillator, as indicated in reference [43]. It should be mentioned that the current work focuses on the intrawell oscillation while more complex dynamic scenarios will be examined in the future work. The size-dependency of the flexoelectric effect on the performance of the energy harvesters in the prebuckling and postbuckling configurations is also investigated. We show the normalized maximum voltage 𝑉̃ = 𝑉max ∕𝑊0 against the beam thickness with or without considering the flexoelectric effect in Fig. 8. The amplitude of the base excitation is the same as that in Fig. 3, i.e. W0 /h = 0.25%. It is seen that the normalized maximum voltage is independent of the beam thickness when only piezoelectric effect is considered, and the normalized value is slightly larger in the postbuckling configuration than that in the prebuckling

International Journal of Mechanical Sciences 173 (2020) 105473

Fig. 8. Variations of the normalized maximum voltage of the energy harvesters with the beam thickness with or without flexoelectric effect.

configuration, which is consistent with the results shown in Fig. 3(a) and (b). However, when the flexoelectric effect is incorporated, the normalized maximum voltage increases significantly with the decrease of the beam thickness from micro to nanoscale, indicating the prominent influence of flexoelectricity on the performance of the energy harvesters at reduced scale. Besides, the normalized maximum voltage due to the combined piezoelectric and flexoelectric effect is larger in the prebuckling configuration than that in the postbuckling configuration when h < 1.5𝜇m. The observed phenomenon is consistent with those shown in Fig. 3(e) and (f), indicating the flexoelectric effect is more prominent than the piezoelectric effect in determining the performance of the energy harvester in the prebuckling configuration. The flexoelectric effect diminishes and both curves (solid line) approach to constant values due to the piezoelectric effect when the beam thickness becomes much larger. 5. Conclusions This work proposes novel nonlinear piezoelectric energy harvesters with flexoelectricity under an axial load and the base excitations. The proposed energy harvester could operate in either prebuckling or postbuckling configuration. Based on the theory of flexoelectricity, the nonlinear electromechanical coupling equations of the proposed energy harvesters under harmonic base excitations are derived utilizing the Hamilton’s principle. The nonlinear Euler–Lagrange discrete equations are obtained by employing the Galerkin’s method and numerically solved by the Runge–Kutta method. Case studies are provided to show the steadystate voltage and power outputs of the energy harvesters in the prebuckling and postbuckling configurations, respectively. It is interesting to find that the frequency response curves for the maximum voltage output of the energy harvesters display a hardening nonlinear behavior in the prebuckling configuration while a softening nonlinear behavior in the postbuckling configuration. The proposed energy harvesters could thus work in a broad bandwidth of frequencies. It is also found that flexoelectricity contributes larger outputs of the energy harvesters in the prebuckling configuration while piezoelectricity induced outputs are more obvious in the postbuckling configuration. Moreover, it is found that load resistance, damping coefficient and the amplitude of the base excitation have significant influences on the performance of the proposed energy harvesters, and both intrawell and interwell oscillations of the energy harvesters could occur in the postbuckling configuration depending on the amplitude of the base excitation. To trigger nonlinear electromechanical responses, a smaller damping coefficient and a larger amplitude of the base excitation are generally required. Specifically, a

Y. Chen and Z. Yan

typical softening nonlinear behavior in the postbuckling configuration is easier to be triggered. Furthermore, the size-dependent property of the flexoelectric effect on the performance of the energy harvesters is demonstrated, which shows the normalized maximum voltage could be prominently enhanced at reduced scale due to flexoelectricity. The results obtained from the current work could be beneficial for designing flexoelectricity-based energy harvesters with optimal performance. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. CRediT authorship contribution statement Yunbin Chen: Methodology, Software, Investigation, Writing - original draft. Zhi Yan: Conceptualization, Investigation, Writing - review & editing, Supervision. Acknowledgments This work is supported by the National Natural Science Foundation of China (No. 11502084). References [1] Mitcheson PD, Miao P, Stark BH, Yeatman EM, Holmes AS, Green TC. MEMS electrostatic micropower generator for low frequency operation. Sensor Actuat A-Phys 2004;115(2–3):523–9. [2] Glynne-Jones P, Tudor MJ, Beeby SP, White NM. An electromagnetic, vibration-powered generator for intelligent sensor systems. Sensor Actuat A-Phys 2004;110(1–3):344–9. [3] Anton SR, Sodano HA. A review of power harvesting using piezoelectric materials (2003–2006). Smart Mater Struct 2007;16(3):R1–R21. [4] Erturk A, Inman DJ. A distributed parameter electromechanical model for cantilevered piezoelectric energy harvesters. J Vib Acoust 2008;130(4):041002. [5] Erturk A, Inman DJ. An experimentally validated bimorph cantilever model for piezoelectric energy harvesting from base excitations. Smart Mater Struct 2009;18(2):025009. [6] Xu R, Lei A, Dahl-Petersen C, Hansen K, Guizzetti M, Birkelund K, Thomsen EV, Hansen O. Fabrication and characterization of MEMS-based PZT/PZT bimorph thick film vibration energy harvesters. J Micromech Microeng 2012;22(9):094007. [7] Kuo CL, Lin SC, Wu WJ. Fabrication and performance evaluation of a metal-based bimorph piezoelectric MEMS generator for vibration energy harvesting. Smart Mater Struct 2016;25(10):105016. [8] Wang ZL, Song J. Piezoelectric nanogenerators based on zinc oxide nanowire arrays. Science 2006;312(5771):242–6. [9] Briscoe J, Dunn S. Piezoelectric nanogenerators – a review of nanostructured piezoelectric energy harvesters. Nano Energy 2015;14:15–29. [10] Cottone F, Vocca H, Gammaitoni L. Nonlinear energy harvesting. Phys Rev Lett 2009;102(8):080601. [11] Erturk A, Hoffmann J, Inman DJ. A piezomagnetoelastic structure for broadband vibration energy harvesting. Appl Phys Lett 2009;94(25):254102. [12] Tang LH, Yang YW. A nonlinear piezoelectric energy harvester with magnetic oscillator. Appl Phys Lett 2012;101(9):094102. [13] Kim P, Seok J. A multi-stable energy harvester: dynamic modeling and bifurcation analysis. J Sound Vib 2014;333(21):5525–47. [14] Zhou S, Zuo L. Nonlinear dynamic analysis of asymmetric tristable energy harvesters for enhanced energy harvesting. Commun Nonlinear Sci Numer Simulat 2018;61:271–84.

International Journal of Mechanical Sciences 173 (2020) 105473 [15] Jiang JX, Liu SG, Zhao D, Feng LF. Broadband power generation of piezoelectric vibration energy harvester with magnetic coupling. J Intel Mat Syst Str 2019;30(15):2272–82. [16] Masana R, Daqaq MF. Electromechanical modeling and nonlinear analysis of axially loaded energy harvesters. J Vib Acoust 2011;133(1):011007. [17] Masana R, Daqaq MF. Relative performance of a vibratory energy harvester in monoand bi-stable potentials. J Sound Vib 2011;330(24):6036–52. [18] Cottone F, Gammaitoni L, Vocca H, Ferrari M, Ferrari V. Piezoelectric buckled beams for random vibration energy harvesting. Smart Mater Struct 2012;21(3):035021. [19] Cottone F, Mattarelli M, Vocca H, Gammaitoni L. Effect of boundary conditions on piezoelectric buckled beams for vibrational noise harvesting. Eur Phys J-Spec Top 2015;224(14–15):2855–66. [20] Zhu Y, Zu J, Su W. Broadband energy harvesting through a piezoelectric beam subjected to dynamic compressive loading. Smart Mater Struct 2013;22(4):045007. [21] Zhu Y, Zu JA. A magnet-induced buckled-beam piezoelectric generator for wideband vibration-based energy harvesting. J Intel Mat Syst Str 2014;25(14):1890–901. [22] Panyam M, Daqaq MF. A comparative performance analysis of electrically optimized nonlinear energy harvesters. J Intel Mat Syst Str 2016;27(4):537–48. [23] Li HT, Qin WY, Zu J, Yang ZB. Modeling and experimental validation of a buckled compressive-mode piezoelectric energy harvester. Nonlinear Dynam 2018;92(4):1761–80. [24] Yang K, Fei F, An H. Investigation of coupled lever-bistable nonlinear energy harvesters for enhancement of inter-well dynamic response. Nonlinear Dynam 2019;96(4):2369–92. [25] Tran N, Ghayesh MH, Arjomandi M. Ambient vibration energy harvesters: a review on nonlinear techniques for performance enhancement. Int J Eng Sci 2018;127:162–85. [26] Yudin PV, Tagantsev AK. Fundamentals of flexoelectricity in solids. Nanotechnology 2013;24(43):432001. [27] Deng Q, Kammoun M, Erturk A, Sharma P. Nanoscale flexoelectric energy harvesting. Int J Solids Struct 2014;51(18):3218–25. [28] Wang KF, Wang BL. An analytical model for nanoscale unimorph piezoelectric energy harvesters with flexoelectric effect. Compos Struct 2016;153:253–61. [29] Yan Z. Modeling of a nanoscale flexoelectric energy harvester with surface effects. Physica E 2017;88:125–32. [30] Basutkar R. Analytical modelling of a nanoscale series-connected bimorph piezoelectric energy harvester incorporating the flexoelectric effect. Int J Eng Sci 2019;139:42–61. [31] Wang KF, Wang BL. Non-linear flexoelectricity in energy harvesting. Int J Eng Sci 2017;116:88–103. [32] Dai HL, Yan Z, Wang L. Nonlinear analysis of flexoelectric energy harvesters under force excitations. Int J Mech Mater Des 2019:1–15. https://doi.org/10.1007/s10999-019-09446-0. [33] Kumar A, Sharma A, Vaish R, Kumar R, Jain SC. A numerical study on flexoelectric bistable energy harvester. Appl Phys a-Mater 2018;124(7):483. [34] Rojas EF, Faroughi S, Abdelkefi A, Park YH. Nonlinear size dependent modeling and performance analysis of flexoelectric energy harvesters. Microsyst Technol 2019;25(10):3899–921. [35] Chen LH, Pan SQ, Fei YY, Zhang W, Yang FH. Theoretical study of micro/nano-scale bistable plate for flexoelectric energy harvesting. Appl Phys A-Mater 2019;125(4):242. [36] Reddy JN. Mechanics of laminated composite plates and shells: theory and analysis. CRC press; 2003. [37] Yan Z, Jiang LY. Flexoelectric effect on the electroelastic responses of bending piezoelectric nanobeams. J Appl Phys 2013;113(19):194102. [38] Hu SL, Shen SP. Electric field gradient theory with surface effect for nano-dielectrics. CMC-Comput Mater Con 2009;13(1):63–87. [39] Erturk A, Inman DJ. Piezoelectric energy harvesting. John Wiley & Sons; 2011. [40] Nayfeh AH, Emam SA. Exact solution and stability of postbuckling configurations of beams. Nonlinear Dyn 2008;54(4):395–408. [41] Emam SA. A theoretical and experimental study of nonlinear dynamics of buckled beams PhD thesis. Virginia Tech; 2002. [42] Liang X, Hu S, Shen SP. Effects of surface and flexoelectricity on a piezoelectric nanobeam. Smart Mater Struct 2014;23(3):035020. [43] Harne RL, Wang KW. A review of the recent research on vibration energy harvesting via bistable systems. Smart Mater Struct 2013;22(2):023001.