Zhuangzhuang He

Gongye Zhang
email: gyzhang@seu.edu.cn

Xin Chen

Yu Cong

Shuitao Gu

Jun Hong
email: junhong@seu.edu.cn



Elastic wave harvesting in piezoelectric-defectintroduced phononic crystal microplates

Keywords: Size effect, Piezoelectric energy harvesting, Modified couple stress theory, Defect band, Phononic crystal, Bandgap

The localization of elastic wave at defect in phononic crystals (PnCs) has been applied to design piezoelectric energy harvesting (PEH) devices. However, the earlier researches were based on classical linear elasticity theory, which failed to predict the structural behavior at the micron scale. Based on a modified couple stress theory (MCST) and a mixed finite element method, the physics builder in COMSOL® Multiphysics is used to develop a novel interface that can capture the microstructure-dependent size effect and be applied in the simulation of a PnC-based PEH device. Size dependence of the newly developed model is demonstrated by three sets of models at different sizes. Numerical results show that when the size effect is considered in the model, the frequencies of the bandgap and defect bands increase with decreasing model size compared to the classical theory. Furthermore, the size reduction changes the internal stiffness ratio in the model, which affects the displacement amplitude, output voltage and output electrical power of the PEH device. Therefore, size effects are inevitable in micron-level models. These discoveries would help the design of small-scale PnC-based PEH devices for enhanced energy collection.

Introduction

The phononic crystal (PnC) defect results in a high degree of elastic wave localization and provides an efficient way of energy collection in piezoelectric energy harvesting (PEH) devices, which has attracted much attention [START_REF] Geng | Thermally-induced transitions of multifrequency defect wave localization and energy harvesting of phononic crystal plate[END_REF][START_REF] Hu | Acoustic-elastic metamaterials and phononic crystals for energy harvesting: a review[END_REF][START_REF] Jo | Double defects-induced elastic wave coupling and energy localization in a phononic crystal[END_REF][START_REF] Jo | Revealing defect-mode-enabled energy localization mechanisms of a one-dimensional phononic crystal[END_REF][START_REF] Wang | Acoustic energy harvesting by piezoelectric curved beams in the cavity of a sonic crystal[END_REF][START_REF] Park | Two-dimensional octagonal phononic crystals for highly dense piezoelectric energy harvesting[END_REF]. Bandgap is an extraordinary property of PnCs, and elastic waves in the bandgap frequency range decay rapidly and cannot propagate through the PnCs [START_REF] Zhang | Magnetically tunable bandgaps in phononic crystal nanobeams incorporating microstructure and flexoelectric effects[END_REF][START_REF] Wang | Tunable and Active Phononic Crystals and Metamaterials[END_REF][START_REF] Jin | Investigation of size effect on band structure of 2D nano-scale phononic crystal based on nonlocal strain gradient theory[END_REF][START_REF] Hong | Tunable bandgaps in phononic crystal microbeams based on microstructure, piezo and temperature effects[END_REF][START_REF] Chen | Low-frequency tunable topological interface states in soft phononic crystal cylinders[END_REF]. Furthermore, by changing a single unit cell structure to destroy the periodicity of the PnCs (i.e., a defect), flat passbands (i.e., defect bands) usually appear in the bandgap [START_REF] Sigalas | Defect states of acoustic waves in a two-dimensional lattice of solid cylinders[END_REF][START_REF] Wu | Point defect states in two-dimensional phononic crystals[END_REF][START_REF] Yao | Propagation of bending waves in phononic crystal thin plates with a point defect[END_REF]. The evanescent waves corresponding to the defect band frequencies can then be localized and amplified inside the defect. Due to the wave enhancement phenomenon, the output power of PEH devices placed at the defect can be several times higher compared to the conventional cantilever-type PEH devices [START_REF] Ma | Acoustic energy harvesting enhanced by locally resonant metamaterials[END_REF].

This innovative idea has been extensively researched. Experiments [START_REF] Wu | Acoustic energy harvesting using resonant cavity of a sonic crystal[END_REF] have shown that the power harvesting efficiency for acoustic energy is 625 times better than that without PnCs, while external loading has an effect on PEH devices. The explicit analytical solution of the 1D defect model [START_REF] Jo | Revealing defect-mode-enabled energy localization mechanisms of a one-dimensional phononic crystal[END_REF][START_REF] Jo | An improved analytical model that considers lateral effects of a phononic crystal with a piezoelectric defect for elastic wave energy harvesting[END_REF][START_REF] Jo | An analytical model of a phononic crystal with a piezoelectric defect for energy harvesting using an electroelastically coupled transfer matrix[END_REF] revealed that the bandgap is a necessary condition for the defect band.

Proper geometric parameters (e.g., piezo thickness, piezo diameter and substrate thickness) can enhance energy transfer and conversion [START_REF] Lee | Enhanced energy transfer and conversion for high performance phononic crystal-assisted elastic wave energy harvesting[END_REF]. For ambient vibrational energies in low-frequency ranges [START_REF] Lin | Modeling and field testing of an electromagnetic energy harvester for rail tracks with anchorless mounting[END_REF][START_REF] Deng | A MEMS based piezoelectric vibration energy harvester for fault monitoring system[END_REF], defective local-resonance type PnCs (often referred to as metamaterials) seem to be more suitable relative to Bragg type [START_REF] Cao | Defect-mode-induced energy localization/harvesting of a locally resonant phononic crystal plate: Analysis of line defects[END_REF]. Single defects have a narrow operating bandwidth, while multiple defects (e.g., decoupled double defects [START_REF] Jo | A graded phononic crystal with decoupled double defects for broadband energy localization[END_REF][START_REF] Jo | A phononic crystal with differently configured double defects for broadband elastic wave energy localization and harvesting[END_REF], coupled double defects [START_REF] Jo | Elastic wave localization and harvesting using double defect modes of a phononic crystal[END_REF] and L-shape defects [START_REF] Jo | L-shape triple defects in a phononic crystal for broadband piezoelectric energy harvesting[END_REF]) can broaden the energy harvesting bandwidth. Multiplefield techniques (e.g., thermal tuning [START_REF] Geng | Thermally-induced transitions of multifrequency defect wave localization and energy harvesting of phononic crystal plate[END_REF][START_REF] Geng | Flexural wave manipulation and energy harvesting characteristics of a defect phononic crystal beam with thermal effects[END_REF], electrical tuning [START_REF] Shakeri | Designing a tunable acoustic resonator based on defect modes, stimulated by selectively biased PZT rods in a 2D phononic crystal[END_REF] and magnetic tuning [START_REF] Deng | A Magnetic-Dependent Vibration Energy Harvester Based on the Tunable Point Defect in 2D Magneto-Elastic Phononic Crystals[END_REF]) allow efficient dynamic adjustment of the energy harvesting operating frequency for adaptation to random environmental frequencies. However, most of them focus on the macroscale level.

Piezoelectric energy harvesting and transfer usually serve some low-power devices (e.g., embedded sensors, medical implants [START_REF] Zaid | Contactless energy transfer using acoustic approach -A review[END_REF][START_REF] Chen | Micro triboelectric ultrasonic device for acoustic energy transfer and signal communication[END_REF] and monitoring devices), which requires miniaturization of the model [START_REF] Deng | A MEMS based piezoelectric vibration energy harvester for fault monitoring system[END_REF]. With device miniaturization, numerical simulations and predictions based on classical linear elasticity theory will encounter a great challenge at the micron scale.

At the microscopic scale, the characteristic scale of the internal microstructure (e.g., inclusions, microcracks and lattices) in the material is close to the material geometry, which makes it impossible to ignore the size effect generated by these microstructures on the mechanical properties of the material. Many experiments (including static bending [START_REF] Choi | Micro-cantilever bending tests for understanding size effect in gradient elasticity[END_REF][START_REF] Lam | Experiments and theory in strain gradient elasticity[END_REF][START_REF] Liebold | Comparison of gradient elasticity models for the bending of micromaterials[END_REF][START_REF] Patel | Novel nickle foil micro-bend tests and the need for a relook at length scale parameter's numerical value[END_REF] and vibration [START_REF] Lei | Size-dependent vibration of nickel cantilever microbeams: Experiment and gradient elasticity[END_REF][START_REF] Li | Experimental investigation on sizedependent higher-mode vibration of cantilever microbeams[END_REF]) have shown that thin structures usually exhibit size effects. In order to solve this problem, a series of higher-order theories have been developed to describe this mechanical behavior, such as the couple stress theory [START_REF] Koiter | Couple-Stress in the Theory of Elasticity[END_REF], the strain gradient theory [START_REF] Mindlin | On first strain-gradient theories in linear elasticity[END_REF], the nonlocal elasticity theory [START_REF] Eringen | Nonlocal polar elastic continua[END_REF], the micromorphic elasticity theory [START_REF] Zhang | A non-classical Bernoulli-Euler beam model based on a simplified micromorphic elasticity theory[END_REF] and the surface elasticity theory [START_REF] Gurtin | A continuum theory of elastic material surfaces[END_REF]. In particular, the modified couple stress theory (MCST) [START_REF] Yang | Couple stress based strain gradient theory for elasticity[END_REF], which includes only one additional material constant, has been successfully applied to explain the microstructuredependent size effect in small-scale structures [START_REF] Ma | A microstructure-dependent Timoshenko beam model based on a modified couple stress theory[END_REF][START_REF] Ma | A non-classical Mindlin plate model based on a modified couple stress theory[END_REF][START_REF] Park | Bernoulli-Euler beam model based on a modified couple stress theory[END_REF][START_REF] Reddy | Microstructure-dependent couple stress theories of functionally graded beams[END_REF], including PnCs [START_REF] Zhang | Elastic wave propagation in 3-D periodic composites: Band gaps incorporating microstructure effects[END_REF][START_REF] Zhang | Band gaps for elastic wave propagation in a periodic composite beam structure incorporating microstructure and surface energy effects[END_REF][START_REF] Zhang | Band gaps for wave propagation in 2-D periodic composite structures incorporating microstructure effects[END_REF]. As a result, the bandgap and defect band in PnCs arising from the Bragg scattering will also be affected by the sizes and microstructures at the micron scales. Hence, there is a need to incorporate the modified couple stress to study PnC-based PEH devices at small scales.

The main objective of this work is to systematically study the size effect on PnC-based PEH devices based on the MCST. Proceeding from the Lagrange multiplier method, we use the physics builder in COMSOL® Multiphysics to develop a new interface (called the MCST interface) that can characterize the size effect, which is an extended version of the built-in solid mechanics interface. This interface is then used in conjunction with the preset electrostatics and circuit interfaces to analyze the model. In Section 2, the theoretical part of this interface is given. Section 3 presents the numerical results of three groups of models with different scales to show how model reduction affects energy collection. Conclusions are made in Section 4.

Formulation

The purpose of this section is to introduce the underlying principles of the new interface and to present the technical details of its implementation in COMSOL® Multiphysics. Two parts will be expanded in this section. Firstly, the MCST is reviewed and the weak form of the motion equations is derived based on Hamilton's principle. Then, two new weak equations combining the Lagrange multiplier method are given. They will be used in Section 3 for the band structure calculation and frequency domain analysis of the model.

Modified couple stress theory

Based on the MCST [START_REF] Yang | Couple stress based strain gradient theory for elasticity[END_REF][START_REF] Zhang | Band gaps for wave propagation in 2-D periodic composite structures incorporating microstructure effects[END_REF], the strain energy UT in region Ω can be written as

1 d 2 ( ) ij ij ij T ij U m V        (1)
where σij and mij are, respectively, the Cauchy stress and couple stress tensors, dV is the volume element, and the strain εij and curvature χij tensors are defined by

    , , , , 1 , 2 1
,

2 ij i j j i ij i j j i u u         (2)
with the displacement vector ui, and the rotation vector θi presented by

, 2 , 1 i ijk k j u    ( 3
)
where εijk is the Levi-Civita symbol. Furthermore, the constitutive equations for isotropic materials can be given by [START_REF] Zhang | Elastic wave propagation in 3-D periodic composites: Band gaps incorporating microstructure effects[END_REF] 2 ,

ij kk ij ij       (4a) 2 2 , ij ij m l    (4b)
where λ and μ are the Lamé constants, δij is the Kronecker delta, and l is a material length scale parameter.

The kinetic energy K of the region Ω is

2 1 d , 2 i K u V      (5)
where ρ is the material mass density, and the overhead "•" represents the first derivative to time t.

In addition, the virtual work W is done by the external force as

    d d . i i i i i i i i W f u c V p u q A           (6)
Here, fi and ci are the body force and body couple stress for per unit volume, respectively. pi and qi are the traction and surface couple for per unit area dA acting on ∂Ω.

Based on Hamilton's principle [START_REF] Zhang | A transversely isotropic magneto-electroelastic Timoshenko beam model incorporating microstructure and foundation effects[END_REF], over the time interval [0, T], the weak form of the motion equations can be obtained from

  0 d 0. T T K U W t          (7)
Combining Eqs. [START_REF] Geng | Thermally-induced transitions of multifrequency defect wave localization and energy harvesting of phononic crystal plate[END_REF][START_REF] Wang | Acoustic energy harvesting by piezoelectric curved beams in the cavity of a sonic crystal[END_REF][START_REF] Park | Two-dimensional octagonal phononic crystals for highly dense piezoelectric energy harvesting[END_REF][START_REF] Zhang | Magnetically tunable bandgaps in phononic crystal nanobeams incorporating microstructure and flexoelectric effects[END_REF], with the virtual displacements vanishing at t = 0 and t = T, we get

    d 0, d i i i i ij ij i i i i i i ij ij m u u f u c p u q V A                         ( 8
)
where the overhead "••" represents the second derivative to time t.

Mixed finite elements for modified couple stress theory

Since the weak form in Eq. ( 8) contains the second-order derivatives of the displacements, the displacements ought to be C 1 -continuous. However, this is a well-known issue. To develop a finite element method for the strain gradient theory, Shu et al. [START_REF] Shu | Finite elements for materials with strain gradient effects[END_REF] proposed a mixed finite element approach using two sets of C 0 -continuous dependent variables (i.e., displacement and displacement gradient) for discretization and Lagrange multipliers to limit the kinematic relationship between them. It was extended by Kwon and Lee [START_REF] Kwon | A mixed element based on Lagrange multiplier method for modified couple stress theory[END_REF][START_REF] Kwon | Three dimensional elements with Lagrange multipliers for the modified couple stress theory[END_REF] to the MCST. Here we briefly review the theoretical part of the method to ensure the completeness of this paper, and more deductive details can be obtained from these studies [START_REF] Shu | Finite elements for materials with strain gradient effects[END_REF][START_REF] Kwon | A mixed element based on Lagrange multiplier method for modified couple stress theory[END_REF][START_REF] Kwon | Three dimensional elements with Lagrange multipliers for the modified couple stress theory[END_REF][START_REF] Deng | Mixed Finite Elements for Flexoelectric Solids[END_REF][START_REF] Nanthakumar | Topology optimization of flexoelectric structures[END_REF][START_REF] Choi | A 10-node tetrahedral element with condensed Lagrange multipliers for the modified couple stress theory[END_REF][START_REF] Mao | Mixed finite-element formulations in piezoelectricity and flexoelectricity[END_REF].

The weak form dependent variables in Eq. ( 8) ultimately point to the displacement. But two sets of independent dependent variables, the displacement vector ui and the relaxed rotation vector i, are considered in the mix element, and the constraint between the relaxed rotation vector i and the physical rotation vector θi is implemented by a Lagrange multiplier vector αi:

( ) d 0. i i i V        (9)
Here, the above equation shows that the connection between the relaxed rotation vector and physical rotation vector is satisfied in the volume-averaged sense [START_REF] Shu | Finite elements for materials with strain gradient effects[END_REF].

For harmonic traveling wave analysis, Eq. ( 8) becomes, with the help of Eq. ( 9) and ci

= 0 = fi = pi = qi 2 ( ) d 0, i ij ij ij ij i i i i i i i i m u u V                                 ( 10
)
where ω is the angular frequency. By applying the Bloch periodic boundary conditions [START_REF] Krattiger | Generalized Bloch mode synthesis for accelerated calculation of elastic band structures[END_REF], the relevant bandgap of the PnCs can be obtained through Eq. [START_REF] Hong | Tunable bandgaps in phononic crystal microbeams based on microstructure, piezo and temperature effects[END_REF].

Considering the damping effect, modifying Eq. ( 10) yields

2 ( ( 1
) ) d 0, s i i i ij ij ij ij i i i i i i m u u V i                                  ( 11
)
where i is an imaginary unit, ηs is the isotropic structural loss factor. Eq. ( 11) will be applied to obtain the frequency response curve of the model.

In COMSOL® Multiphysics, we use the physics builder to package the above into an MCST interface. For the subsequent numerical simulations, the displacement vector is discretized using the Lagrange quadratic element, while the relaxed rotation and Lagrange multiplier vectors are discretized using the Lagrange linear element, as shown in Fig. 1. The 3×3 and 2×2 Gaussian quadrature methods are utilized for displacement and relaxed rotation, respectively. In addition, the validation of the MCST interface developed in this section is illustrated in Appendix A with a pure shear example. Fig. 1 A hexahedral element with 27 nodes for MCST interface, with the hollow circle being the displacement vector ui and the solid circle representing the displacement ui, the relaxed rotation i and Lagrange multiplier αi vectors.

PnC-based PEH device

What exactly are the size effects on PnC-based PEH devices? For this purpose, numerical results for three different sets of models at different scales are given in this section. Compared with classical theory, the couple stress effect is revealed by comparing the band structure, defect band and PEH device's performance, respectively. As shown in Fig. 2(a), the PnC unit cell consists of a hard iron core (white part) and a soft epoxy substrate (light green part). The core is centered and penetrates the substrate, while the cell geometry is symmetric about the neutral layer of the substrate. In this section, models are presented for three different body epoxy plate thicknesses (i.e., h = 1.5 mm, 150 μm and 100 μm) to demonstrate the size effect. In addition, the lattice constant of the unit cell a = 20h, the length of the iron core a1 = 10h, and the core thickness tir = 4h. The material properties of the iron and epoxy are given in Table 1.

Bandgap analysis for unit cell

Table 1. Material properties of iron and epoxy [START_REF] Zhang | Microstructure-dependent Band Gaps for Elastic Wave Propagation in a Periodic Microbeam Structure[END_REF]. and for epoxy bh = 24 µm [START_REF] Lam | Experiments and theory in strain gradient elasticity[END_REF].

Bloch periodic boundary conditions are arranged around the unit cell and the reciprocal wave vector is scanned along the boundary line of the irreducible first Brillouin zone (see From Fig. 3, all unit cells have a complete bandgap, which means that all elastic waves (i.e., both in-plane and out-of-plane waves) of frequencies in the bandgap range will be filtered.

When h = 1.5 mm, the current model predicts the bandgap values that are consistent with those of the classical model. However, the bandgap frequency of the current model is higher than that of the classical model when h = 150 µm, and this difference is further reinforced at h = 100 µm. This finding implies that the size effect on the bandgap is significant only in very thin plates [START_REF] Zhang | Band gaps for flexural elastic wave propagation in periodic composite plate structures based on a non-classical Mindlin plate model incorporating microstructure and surface energy effects[END_REF]. In addition, considering the microstructure-dependent size effect stiffens the model and leads to an elevated eigenfrequency, as expected [START_REF] Lei | Size-dependent vibration of nickel cantilever microbeams: Experiment and gradient elasticity[END_REF][START_REF] Li | Experimental investigation on sizedependent higher-mode vibration of cantilever microbeams[END_REF]. 1.

Frequency (Hz) Frequency (Hz)

Frequency (Hz)

Frequency (Hz)

Defect bandgap analysis

The periodicity of PnCs is the main factor for the bandgap generation, and if this periodicity is disrupted by locally transforming the material properties or structure inside PnCs (i.e., introducing a defect), the elastic wave will localize at the defect. The major focus in this subsection is on how the size effect impacts the defect band frequencies. 2. Note that the material properties of the PZT-4D in Table 2 are directly taken from the COMSOL built-in material library.

Furthermore, at the micro-scale, the size effect of ceramic materials is not significant and can usually be ignored (l = 0) [START_REF] Choi | Micro-cantilever bending tests for understanding size effect in gradient elasticity[END_REF]. The circuit of the PEH devices is depicted in Fig. 4(b). Note that the top and bottom surfaces of the two piezoelectric sheets are considered to be completely covered by the electrodes, while the electrode thickness is negligible with respect to the whole model. The top and bottom electrodes of the PEH devices are connected to each other and to the input section of the resistor, and the remaining electrodes are connected to the output section of the resistor and grounded. The external resistor is uniformly set to 100 MΩ for simulating an open-circuit condition.

For the case of PEH devices, the electrostatics and electrical circuit interfaces in COMSOL® Multiphysics are used in conjunction with the solid mechanics interface (i.e., the classical model) or the MCST interface (i.e., the current model) for calculating the defect band frequencies. The electrostatics interface is coupled with the solid mechanics or MCST interface to depict the piezoelectric effect. The electrical circuit interface is responsible for implementing the circuit configuration in Fig. 4(b) and determining the performance of the PEH devices.

Under the open-circuit condition (i.e., resistance of 100 MΩ), Fig. 5 shows the defect bands of the current and classical models within the bandgap at different scales (i.e., h = 1.5 mm, 150µm, 100 µm). In Appendix B, we provide all the defect mode shapes (displacement magnitudes) in full and summarize them into three types (i.e., Modes A, B and C) listed in Fig. 6. Referring to the previous study by Lee et al. [START_REF] Lee | Enhanced energy transfer and conversion for high performance phononic crystal-assisted elastic wave energy harvesting[END_REF][START_REF] Zhang | Band gaps for flexural elastic wave propagation in periodic composite plate structures based on a non-classical Mindlin plate model incorporating microstructure and surface energy effects[END_REF], the target bandgap and defect band frequencies can be obtained by adjusting the mass or stiffness of the unit cell or PEH device.

The information given in Fig. 5 can be summarized as follows: for the same type of defect mode, the current model always predicts a higher frequency than the classical model. As the plate thickness h decreases, this difference becomes larger. In particular, when h = 100 µm, the corresponding frequency of Mode C disappears in the bandgap. It is also important to emphasize that the size effect can be neglected in the large-scale case (with h = 1.5 mm here) since both models predict the same defect band frequencies. From Fig. 6, we confirm that the defect state modes for all models are out-of-plane versions. Mode B, as a unipolar class defect mode shape, is the best deformation for PEH device energy collection efficiency. This is mainly attributed to two reasons: (1) the PEH device should have sufficient deformation to generate electro-elastic coupling for output electrical energy; and (2) the PEH device should avoid being deformed in both tension and compression to avoid voltage cancellation. In conclusion, for all modes, the frequencies corresponding to Mode B will be the focus of subsequent analysis. 1 and2. For all frequency domain analyses, the mechanical damping of the entire structure is set to a structured loss factor ηs = 0.0001. The maximum grid size is always less than 1.67 times the substrate thickness, which is about one-tenth of the corresponding excitation wavelength and is sufficient to achieve highly accurate numerical results. It is important to emphasize that the grid setting strategy is valid for all cases in this paper. We also provide specific meshing schemes in Appendix B to ensure reproducibility in the simulations.

The frequency responses of the z-displacement and output voltage for the current and classical models with three thickness cases are plotted in Fig. 8, and the parameter peaks are recorded in Table 3. The displacement amplitude is the value at the top center of the upper PEH device, and the output voltage comes from the potential difference between the external load input and the output section in the electrical circuit interface. Nearly identical peak displacement (137 nm) and peak output voltage (0.50 V) are obtained for the current and classical models at h = 1.5 mm. When h = 150 and 100 μm, the peak displacement and output voltage of the classical model are guaranteed to be consistent with the part of h = 1.5 mm. On the contrary, the displacement and voltage of the current model at h = 150 μm are raised to 147 μm and 0.56 V. The difference is further amplified at h = 100 μm (with the peak be observed that the resonant frequency of the current model is gradually lifted relative to the classical model, which is consistent with the conclusions of Section 3.2.

It should be noted that the higher peak displacements and voltages obtained for the current model in the small-scale case (with h = 150 and 100 μm here) are due to the change in the stiffness ratio of the PnCs to the PEH devices (i.e., the effective stiffness of the PnCs increases for the current model relative to the classical model, while the PEH one remains constant). In addition, higher output voltages can be obtained by connecting two PEH devices in series [START_REF] Jo | Elastic wave localization and harvesting using double defect modes of a phononic crystal[END_REF].

Optimal external load for maximum electrical output power

For all cases, the external resistance is swept from 0 to 100 MΩ at the resonant frequency (see Table 3), and the rest of the settings are consistent with Section 3.3. Fig. 9 shows the output voltage and output power corresponding to different resistances at the resonant frequency. According to Ohm's law, the electric output power is defined as voltage 2 /resistance.

From Figs. 9(a, c, e), the output voltage increases monotonically and converges to a peak for all models as the external resistance increases. In all cases, the output voltage of the current model is always higher than that of the classical model, regardless of the external resistance. This difference is further amplified as the scale decreases.

Figs. 9(b, d, f) show the peak output power of each model and the corresponding optimal external resistance. The optimal electrical output for all models is summarized in Table 4. It can be shown that, regardless of plate thickness, the classical model has a constant value of the optimal resistance (with 98.1kΩ) and peak output power (0.65μW) under the resonant frequency wave excitation. While for the current model, at h = 1.5 mm, the same value results are obtained with the classical model. For h = 150 μm, the optimal resistance and peak output power of the current model are 100.2 kΩ and 0.8 μW, respectively, which are both larger than the values of the classical model. The numerical difference is further stretched at h = 100 μm with 103.4 kΩ and 0.87 μW, respectively. This indicates that, in contrast to classical continuum theory, the influence of the microstructure-dependent size effect can be very important at small scales, boosting the output voltage and power of the PnC-based PEH device. This finding further supports the miniaturization of the PnC-based PEH device.

Conclusion

In this paper, we introduce the microstructure-dependent size effect in the PnC-based PEH device, which provides a better description of the small-scale model mechanical and electrical behaviors. First, the MCST interface that can characterize the size effect was developed using the Lagrange multiplier method in COMSOL® Multiphysics. Second, using the newly developed interface, the band structure, defect band and PEH device's performance predicted by the current model were compared to that predicted by the classical model for three sets of scale size PnC-based PEH devices. The newly developed model can be applied to predict results for PnC-based PEH devices at all scales, while the classical model can only predict the relevant solutions at large scales. When the device is not very small, the current model predicts the same bandgap frequencies as the classical model.

Considering the size effect raises the bandgap and defect band frequencies in the current model compared to the classical model. Especially, the effect on defect band frequency is significant. Moreover, this difference is further amplified as the model scale decreases. The frequency domain analysis of the open-circuit condition shows that the displacement, voltage, and peak output power of the current model are all elevated with respect to the classical model due to the change in the stiffness ratio between the PnCs and PEH devices. In addition, a unipolar class defect mode shape can always be detected in both current and classical models for the best PnC-based PEH device energy collection efficiency.

When the scale comes to the micron level, models have higher operating frequencies. The predictions and numerical simulations for the PnC-based PEH device by the classical linear elasticity theory have great limitations, and it is necessary to consider the size effect. Furthermore, energy harvesting at a micron or even smaller scale can be realized and enhanced by the current PnC-based PEH device. 1.The analytical solution is adopted from Park et al. [START_REF] Park | Variational formulation of a modified couple stress theory and its application to a simple shear problem[END_REF]. The classical solution corresponding to l = 0 is also presented.

Appendix B

B.1. Meshing scheme

The maximum mesh length for all FE models in this paper is no more than 1.67 times the thickness of the model epoxy substrate. The specific division scheme is shown in

Fig. 2

 2 Fig.2 Schematic diagram of the PnC: (a) unit cell consists of a core with the length a1 and thickness tir, and a substrate with the lattice constant a and thickness h, and (b) irreducible first Brillouin zone.

Fig. 2 (

 2 Fig.2(b)). The band structure of the unit cell is obtained by calculating the eigenfrequencies. In this paper, the model incorporating size effect (i.e., modeled with the MCST interface) is considered as the current model, while the classical model is modeled with the solid mechanics interface built in COMSOL. The displacements in the solid mechanics interface are Lagrange quadratic discretization, while the discretization strategy in the MCST interface we have accounted for in Section 2. The bandgaps of the three scaled unit cells (i.e., h = 1.5 mm, 150 μm and 100 μm) are plotted in Fig 3.

Fig. 3

 3 Fig. 3 Band gaps for the PnC predicted by: (a) the current model with h = 1.5 mm, (b) the classical model with h = 1.5 mm, (c) the current model with h = 150 μm, (d) the classical model with h = 150 μm, (e) the current model with h = 100 μm, and (f) the classical model with h = 100 μm. The lattice constant of the unit cell a = 20h, the length of the iron core a1 = 10h, the core thickness tir = 4h, and the material properties employed in the calculation are tabulated in Table1.

Fig. 4 A

 4 Fig. 4 A PnC-based PEH system: (a) a geometric configuration of a 5×5 defected PnCs plate with two piezoelectric sheets in the middle position, with the radius rpz and thickness tpz, and (b) circuit diagram for two piezoelectric sheets.

Fig. 5

 5 Fig. 5 Defect bands for the PnC-based PEH device predicted by: (a) the current model with h = 1.5 mm, (b) the classical model with h = 1.5 mm, (c) the current model with h = 150 μm, (d) the classical model with h = 150 μm, (e) the current model with h = 100 μm, (f) the classical model with h = 100 μm. The lattice constant of the unit cell a = 20h, the length of the iron core a1 = 10h, the core thickness tir = 4h, the radius of the piezoelectric sheet rpz = 2.5h, thickness tpz = 0.5h, the external resistor is100 MΩ, and the material properties employed in the calculation are tabulated in Tables1 and 2.

Fig. 6

 6 Fig. 6 Three types of defect mode shapes (displacement magnitudes) of the current and classical models at different scales: (a) Mode A, (b) Mode B, and (c) Mode C.

3. 3 Fig. 7

 37 Fig. 7 Schematic diagram of the frequency domain analysis model, with the excitation to be the distance of 5a/3 from the front-most PnCs. Two perfectly matched layers (PMLs) are set at the substrate boundary in the x-direction, and continuous periodic boundary conditions (PBCs) are imposed in the y-direction.

Fig. 8

 8 Fig.8 The PEH performance of the current and classical models: (a) z-displacement at h = 1.5 mm, (b) output voltage at h = 1.5 mm, (c) z-displacement at h = 150 μm, (d) output voltage at h = 150 μm, (e) z-displacement at h = 100 μm, and (f) output voltage at h = 100 μm.

Fig. 9

 9 Fig. 9 Variation of outputs concerning external resistance of the current and classical models: (a) output voltage at h = 1.5 mm, (b) output power at h = 1.5 mm, (c) output voltage at h = 150 μm, (d) output power at h = 150 μm, (e) output voltage at h = 100 μm, and (f) output power at h = 100 μm. The electric output power is defined as voltage 2 /resistance.

Fig. A. 2

 2 Fig. A.2 Comparison of three numerical results (i.e., u1) at the cut line (L/2, y, 0) for the pure shear problem. The FE model is calculated based on the MCST interface developed in Section 2.2. The length L = 300 μm, width b = 150 μm, thickness h = 30 μm, and the material properties of the epoxy employed in the calculation are tabulated in Table 1.The analytical

FigFig. B. 1

 1 Fig. B.1 Meshing schemes for all models mentioned in this paper: (a) unit cell, (b) supercell, and (c) frequency domain model.

Fig. B. 2

 2 Fig. B.2 Displacement defect mode shapes for the current model on (a) 8.76 kHz at h =1.5 mm, (b) 9.59 kHz at h =1.5 mm, (c) 10.40 kHz at h =1.5 mm, (d) 89.79 kHz at h =150 μm, (e) 98.44 kHz at h =150 μm, (f) 106.80 kHz at h =150 μm, (g) 138.73 kHz at h =100 μm, and (h) 152.25 kHz at h =100 μm.

Table 2 .

 2 Material properties of PZT-4D.

	Material properties	ρ(kg/m 3 )	c11(GPa)	c12(GPa)	c13(GPa)	c33(GPa)	c44(GPa)	e31(C/m 2 ) e33(C/m 2 ) e15(C/m 2 )	ϵ11(10 -9 F/m)	ϵ33(10 -9 F/m)
	PZT-4D	7600	153.83	98.46	93.10	128.24	23.81	-4.7303	15.2586	13.0952	7.05	6.75

Table 3 .

 3 Peak z-displacement and peak voltage for the current and classical models at three different scales (h = 1.5 mm, 150 μm and 100 μm).

			Resonance	Peak displacement	Peak voltage
	Plate thickness h	Model
			frequency (kHz)	(nm)	(V)
	1.5 mm	Classical	9.59	137	0.50
	1.5 mm	Current	9.59	137	0.50
	150 μm	Classical	95.89	137	0.50
	150 μm	Current	98.44	147	0.56
	100 μm	Classical	143.84	137	0.50
	100 μm	Current	152.24	149	0.60

Table 4 .

 4 Optimal external resistance and peak output power for the current and classical models at three different scales (h = 1.5 mm, 150 μm and 100 μm).

	Voltage (V)	Output Power (μW)
		Optimal external	Peak output power
	Plate thickness h	Model
		resistance (kΩ)	(μW)

Frequency (Hz)

Acknowledgements

This work was supported by the National Natural Science Foundation of China (Grant No. 12002086), the Zhishan Youth Scholar Program of SEU, and the Fundamental Research Funds for the Central Universities.

Appendix A

It is shown here that the MCST interface developed in Section 2.2 is validated with a pure shear example. Park et al. [START_REF] Park | Variational formulation of a modified couple stress theory and its application to a simple shear problem[END_REF] provide an analytical solution to a pure shear problem for the MCST. With the help of the MCST interface, we compare the numerical values of the finite element (FE) model to them.

As shown in

The remaining parts are the traction free boundary conditions. It should be noted that the premise of the analytical solution proposed by Park et al. [START_REF] Park | Variational formulation of a modified couple stress theory and its application to a simple shear problem[END_REF] assumes that L, b→∞, which is the reason for imposing symmetry boundary conditions (i.e., constraints on the z = 0, b planes) in the FE model.