On the E ﬀ ects of Structural Coupling on Piezoelectric Energy Harvesting Systems Subject to Random Base Excitation

: The current research investigates the novel approach of coupling separate energy harvesters in order to scavenge more power from a stochastic point of view. To this end, a multi-body system composed of two cantilever harvesters with two identical piezoelectric patches is considered. The beams are interconnected through a linear spring. Assuming a stochastic band limited white noise excitation of the base, the statistical properties of the mechanical response and those of the generated voltages are derived in closed form. Moreover, analytical models are derived for the expected value of the total harvested energy. In order to maximize the expected generated power, an optimization is performed to determine the optimum physical and geometrical characteristics of the system. It is observed that by properly tuning the harvester parameters, the energy harvesting performance of the structure is remarkably improved. Furthermore, using an optimized energy harvester model, this study shows that the coupling of the beams negatively a ﬀ ects the scavenged power, contrary to the e ﬀ ect previously demonstrated for harvesters under harmonic excitation. The qualitative and quantitative knowledge resulting from this analysis can be e ﬀ ectively employed for the realistic design and modelling of coupled multi-body structures under stochastic excitations. H.M. (Hamid Moeenfard) and H.H.K.; methodology, H.M. (Hamidreza Masoumi), H.M. (Hamid Moeenfard) and H.H.K.; software, H.M. (Hamidreza Masoumi).; validation, H.M. (Hamidreza Masoumi); formal analysis, H.M. (Hamidreza Masoumi), H.M. (Hamid Moeenfard) and H.H.K.; investigation, H.M. (Hamidreza Masoumi), H.M. (Hamid Moeenfard) and H.H.K.; resources, H.M. (Hamidreza Masoumi), H.M. (Hamid Moeenfard), M.I.F. and H.H.K.; data curation, H.M. (Hamidreza Masoumi), H.M. Moeenfard) and H.H.K.; writing—original H.M. Masoumi), H.M. Moeenfard); and H.H.K. and M.I.F.; visualization, H.M. (Hamidreza Masoumi), H.M. Moeenfard) and H.H.K.; H.M. Moeenfard) and H.H.K. project administration, H.M. Moeenfard); Moeenfard) H.H.K. authors


Introduction
Extracting electrical energy from ambient vibration has received extensive attention in recent years. The use of this source of energy can potentially eliminate the need for battery replacement [1]. The use of thermoelectric [2], electromagnetic [3], electrostatic [4], and piezoelectric [5] effects are the most significant for energy harvesting, and vehicle and machine vibrations are some of the most common energy sources that can be employed for scavenging energy [6]. The power generated from a sample excitation is approximately proportional to the cube of the excitation frequency [7]. Hence, extracting sufficient electrical power from low frequency excited media is a major challenge. By matching the resonance frequencies of the system with that of the excitation, many researchers have devised appropriate techniques to tackle this limitation. As an example, Challa et al. [8] presented the idea of a tunable resonance frequency using magnetic force which is capable of broadening the resonance frequency by ±20% its untuned value. Tehrani and Elliot [9] utilized a cubic nonlinear

Mathematical Modelling of the System
In this section, the differential equations governing the dynamic behavior of the system are derived using Hamilton's principle. These equations and the corresponding boundary conditions are then further employed to find the natural frequencies and mode shapes of the system. The schematic view of the energy harvester under study, consisting of two different beams connected via a linear spring, is shown in Figure 1. The beams are clamped at one end while connected through a linear spring at the other. Two similar piezoelectric patches are attached to each beam while the whole system is under a random base excitation.
Aerospace 2020, 7, x FOR PEER REVIEW 3 of 28 In this section, the differential equations governing the dynamic behavior of the system are derived using Hamilton's principle. These equations and the corresponding boundary conditions are then further employed to find the natural frequencies and mode shapes of the system. The schematic view of the energy harvester under study, consisting of two different beams connected via a linear spring, is shown in Figure 1. The beams are clamped at one end while connected through a linear spring at the other. Two similar piezoelectric patches are attached to each beam while the whole system is under a random base excitation.  Figure 1 shows that the lengths of both piezoelectric patches are very small compared to the length of the beams. As a result, they have negligible effects on the mode shapes. Therefore, in the upcoming derivations (which aim to find the mode shapes of the coupled system) the effects of the piezoelectric patches are ignored.

Eigen-Value Problem
By utilizing Hamilton's principal, the normalized free vibration equations of the system and the corresponding kinematic and natural boundary conditions are derived as [21] Equations of motion:   Figure 1 shows that the lengths of both piezoelectric patches are very small compared to the length of the beams. As a result, they have negligible effects on the mode shapes. Therefore, in the upcoming derivations (which aim to find the mode shapes of the coupled system) the effects of the piezoelectric patches are ignored.

Eigen-Value Problem
By utilizing Hamilton's principal, the normalized free vibration equations of the system and the corresponding kinematic and natural boundary conditions are derived as [21].
Equations of motion: Kinematic boundary conditions: Natural boundary conditions: (1,t) = 0 (4) (1,t) = 0, i = 1, 2 In these equations L s , is the length of the beams, E s is the Young's modulus of the beams material, s is the area cross section of the ith beam and k is the elastic constant of the connecting linear spring. Moreover, the normalized variables x, w i and t are defined by whereŵ i is the deflection of the ith beam and T = L 2 s ρ s A (1) s /E s I (1) s in which ρ s is the material density. The dynamic response of a linear continuous system can be expressed as a linear combination of its mode shapes [21]. Therefore, an essential first step in the dynamic modelling of a continuous structure, either single or multi-body, is to find the mode shapes of that system. A mode shape of a system is a state of motion of that system in which all of its elements move with the same frequency and phase, but not necessarily with the same amplitude. Appendix A presents the derivation of the exact natural frequencies and the mode shapes of the system. To verify the accuracy of the proposed analytical technique, a system with physical and geometrical parameters given in Table 1 is considered. Using the numerical values reported in Table 1, the first three normalized natural frequencies of the system are obtained as given in Table 2. In this table, the analytical results are compared with those of finite element (FE) simulations using the commercial software Abaqus. The numerical and analytical results are in excellent agreement. In order to investigate the accuracy of the analytical mode shapes, Figures 2-4 are provided. In these figures, the first three analytical and FE modes are compared, and acceptable agreement is observed. To make the numerical and analytic modes comparable, they have both been normalized such that

Dynamic Analysis
Lagrange's equations are employed here to derive the differential equations governing the dynamic behavior of the system. This requires the analytical expressions for the potential and kinetic energies, as well as the virtual work done on the system. The potential energy of the system shown in Figure 1 is composed of three parts: (1) strain energy of the substructure, , (2) potential energy of the piezoelectric material, , and (3) strain energy of the spring, . • Strain Energy of the Substructure Assuming linear elastic material for the beams, the strain energy of the harvester can be written as where L is the length of the piezoelectric patches (see Appendix B) and • Potential Energy of the Piezoelectric Layers The potential energy of the piezoelectric material, , is composed of a strain energy and a potential electrical component . The strain energy of the piezoelectric layers in the system shown in Figure 1 can be expressed as [22]

Dynamic Analysis
Lagrange's equations are employed here to derive the differential equations governing the dynamic behavior of the system. This requires the analytical expressions for the potential and kinetic energies, as well as the virtual work done on the system. The potential energy of the system shown in Figure 1 is composed of three parts: (1) strain energy of the substructure, π s , (2) potential energy of the piezoelectric material, π p , and (3) strain energy of the spring, π e .
• Strain Energy of the Substructure Assuming linear elastic material for the beams, the strain energy of the harvester can be written aŝ where L p the length of the piezoelectric patches (see Appendix B) and • Potential Energy of the Piezoelectric Layers The potential energy of the piezoelectric material,π p , is composed of a strain energyÛ p and a potential electrical componentŴ ie . The strain energy of the piezoelectric layers in the system shown in Figure 1 can be expressed as [22] Aerospace 2020, 7, 93 7 of 26 where σ (1) p and σ (2) p are respectively the normal stress components of the piezoelectric layers attached to the thick and thin beams, while ε (1) p and ε (2) p are their strain counterparts. Moreover, ∀ (k) p is the volume of the kth piezoelectric layer. The following constitutive equations show the relation between the stress and strain in the piezoelectric layers [23]: In Equations (11) and (12), D   [24] (wherev k (t) is the generated voltage in the kth piezo layer) and employing the stress-strain relationships, one may simplify theÛ p expression asÛ where I The internal electrical energy generated in the piezoelectric layers,Ŵ ie , can be expressed as [23] By substituting D in which Thus, the total potential energy of piezoelectric layers is obtained aŝ Aerospace 2020, 7, 93 8 of 26 • Strain Energy of the Connecting Spring The elastic strain energy of the linear connecting spring,π e , is simplŷ π e = 1 2 k(ŵ 1 (L s ,t) −ŵ 2 (L s ,t)) 2 (21) Assuming that the mass of the spring is negligible, the kinetic energy of the system can be derived by adding the kinetic energies of the constitutive beams and the piezoelectric layers aŝ whereŶ is the base displacement.
To account for the modal damping, it is assumed that the system is vibrating in a linearly damped viscous media. The damping coefficient for both beams is assumed to be c per unit length. In such a condition, the virtual work done on the system due to this damping will be as In addition, as the electrical charge is Q i , i = 1, 2, in the load circuit with resistance R (i) l , the virtual work is δŴ Hence, the effective virtual work done to the system can be summarized as The dynamic response of linear systems can be effectively approximated as a linear combination of their mode shapes. This theory has been examined and verified in many published studies [5,25]. Hence, for the system under study, whereq i represents the contribution of the ith mode in the dynamic response. Now,q i , as well asv 1 and v 2 , can be considered as the generalized coordinates. By substituting Equation (26) into (7), (20), (21), (22), and (25), while employing Lagrange's equation, after some simple manipulations, the governing equations are derived. For notational simplicity, the following normalized variable is defined pc is the distance between the center of the piezoelectric layer and the neutral axis and d 31 and ε s 33 are piezoelectric constants whose typical values are given in Table 3. Using the normalized variable q =q/L s , the dimensionless governing equations are reduced to . In Equation (29), ζ and ϑ k are defined as In these equations, n×1 are functions of the mode shapes as well as the physical and geometrical parameters of the system, and are given by Aerospace 2020, 7, 93 10 of 26 where γ (k) in Equation (35) is calculated as In addition, I (k) p in Equation (36) is In the next step, the dynamic Equations (28) and (29) are transformed into state space. To do so, the state vector → P is defined as where the superscript T denotes the transpose operator. The state space equations governing the dynamic behavior of the system are then ..
To investigate the number of modes that contribute significantly to the dynamic response, a physical system with the specifications given in Table 1 (for the substructure) and Table 3 (for the piezoelectric materials) is considered. Moreover, here and in the rest of this paper, unless otherwise stated, the damping c is selected such that the modal damping of the first mode becomes a nominal value of 0.01.

Validation
To further investigate the accuracy of the presented formulation, the natural frequencies of the harvester are compared with those presented in [26]. To this end, a harvester with zero coupling and mechanical properties as stated in Tables 4 and 5 is considered. The frequency response curves for the sample harvester, resulting from the base acceleration for R l = 1 MΩ (see Table 6) and R l = 100 Ω (see Table 7), are illustrated in Figure 5. The exact values of the natural frequencies in which the voltage response peaks, are compared with the results presented in [26]. Regarding the peak and error values presented in the tables, it is evident that the present model acceptably predicts the resulting voltage from the harvester.  Figure 5. The exact values of the natural frequencies in which the voltage response peaks, are compared with the results presented in [26]. Regarding the peak and error values presented in the tables, it is evident that the present model acceptably predicts the resulting voltage from the harvester.

Random Vibration Analysis
A necessary first step in the random vibration analysis of a dynamic system is to obtain the frequency response functions. The frequency response matrix H ( ..

Y)
P (ω) . of the linear system (Equation (39)), can be derived as The mathematical details to obtain Equation (42) are given in Appendix C and j * = √ −1. The frequency response function corresponding to w k (x, t) can also be easily derived by substituting

Y)
P i (ω) being the ith element of the frequency response matrix derived in Equation (42)) into the normalized form of Equation (26) and cancelling exp( j * ωt) from both sides. Consequently, one can get (43) Figure 6 shows the magnitude of the frequency response functions of the normalized tip displacements w 1 (1, t) and w 2 (1, t) of the undamped system. The frequency responses of the generated voltages are also illustrated in Figure 7. These frequency responses will be used later to derive the statistical properties of the harvested power.
Aerospace 2020, 7, x FOR PEER REVIEW 13 of 28 being the ith element of the frequency response matrix derived in Equation (42)) into the normalized form of Equation (26) and cancelling exp(j * ωt) from both sides. Consequently, one can get Figure 6 shows the magnitude of the frequency response functions of the normalized tip displacements (1, ) and (1, ) of the undamped system. The frequency responses of the generated voltages are also illustrated in Figure 7. These frequency responses will be used later to derive the statistical properties of the harvested power.
Using Laplace transform, this impulse response function can be easily derived as The following equation characterizes the relationship of the frequency and impulse response functions [15] Having established the impulse response function, the dynamic response of the system to a general input .. Y(t) can be obtained using the following convolution integral [15].
Assuming a stationary process for .. (47) can be further simplified as

Y(t), Equation
Y(t). By choosing i = j in Equation (48), the autocorrelation function R p i (τ) is obtained as Having a closed form relation for R p i p j (x, τ), the autocorrelation function for w k (x, t) can also be derived as Upon simplification, R w k (x, τ) can be expressed as The cross spectral density of p i and p j is defined as s p i p j (ω) = 1/2π ∞ −∞ R p i p j (τ) exp(−jωτ)dτ. By substituting Equation (48) into this equation, one gets By defining γ = τ − θ 2 + θ 1 and changing the order of integration, it can be proved that where S .. Y (ω) is the base acceleration spectral density. By simplifying Equation (53) while considering Equation (45), S p i p j (ω) is derived as where H * p i (ω) is the complex conjugate of H p i (ω). By selecting i = j, Equation (54) can be used to derive a compact formula for S p i (ω) as Considering Equation (55) and choosing i = 2n + 1 and i = 2n + 2, S v 1 (ω) and S v 2 (ω) are derived as The spectral density of p i , can be employed to determine the spectral density of w i (x, t). In fact, S w k (x, ω) can be obtained by taking the Fourier transform of R w k (x, τ).
By substituting Equation (51) into (58) and considering the definition of the cross spectral density given in Equation (54), S w k (x, ω) can be simplified as Considering thatσ (k) s = −E s y∂ 2ŵ /∂x 2 , with the same procedure which was used for deriving Equation (59), the following equation can be obtained for the spectral density function of normalized normal stress components of the beams One of the most important parameters which influences the accuracy of the dynamic model is the number of modes included in the model. The sufficient number of modes depends on the type and the frequency extent of the random loads. Here in this study, the excitation is assumed to be band limited white noise with the spectral density shown in Figure 8. The S 0 in this figure is assumed to be 0.0169 (corresponding to a nominal value of 15.5031 m 2 /s 3 ). Moreover, in this figure, a 0 is assumed to be 36.966 (corresponding to nominal value of 50 Hz) which is slightly lower than the fourth natural frequency of the system, 53.678 Hz. Since the maximum excitation frequency is between the third and the fourth natural frequencies of the system, it would be reasonable to include four modes in dynamic analysis. Thus, hereafter in this paper, unless otherwise mentioned, we will consider four modes in the simulations.
Considering that One of the most important parameters which influences the accuracy of the dynamic model is the number of modes included in the model. The sufficient number of modes depends on the type and the frequency extent of the random loads. Here in this study, the excitation is assumed to be band limited white noise with the spectral density shown in Figure 8. The in this figure is assumed to be .0169 (corresponding to a nominal value of 15.5031 m 2 /s 3 ). Moreover, in this figure, is assumed to be 36.966 (corresponding to nominal value of 50 Hz) which is slightly lower than the fourth natural frequency of the system, 53.678 Hz. Since the maximum excitation frequency is between the third and the fourth natural frequencies of the system, it would be reasonable to include four modes in dynamic analysis. Thus, hereafter in this paper, unless otherwise mentioned, we will consider four modes in the simulations. The spectral densities corresponding to the resulting voltage (i.e.,

( ) and
( )) are shown in Figure 9. As the figure shows, both spectral density functions have two peaks, the second of which is larger than the first one. Furthermore, the peaks of ( ) and ( ) occur at the same frequency which is due to the coupled dynamics of beams. The spectral densities will be used later to find analytical expressions for the harvested electrical powers. Furthermore, the spectral The spectral densities corresponding to the resulting voltage (i.e., S v 1 (ω) and S v 2 (ω)) are shown in Figure 9. As the figure shows, both spectral density functions have two peaks, the second of which is larger than the first one. Furthermore, the peaks of S v 1 (ω) and S v 2 (ω) occur at the same frequency which is due to the coupled dynamics of beams. The spectral densities will be used later to find analytical expressions for the harvested electrical powers. Furthermore, the spectral densities S It has to be clarified that since expected maximum stress is zero (i.e., ( ) = 0), it cannot be utilized as a proper criterion to investigate the maximum deflection. Instead, since E σ (t) is of the order of magnitude of ( ), one can use it to provide an acceptable approximation of the order of magnitude of the maximum stress in the beams. So, the mean square stresses given in Equations (63) and (64) can be employed as a criterion to determine the physical failure limits during the random excitation of the system.

Parametric Study
In order to investigate the effects of different design parameters on the harvested electrical power, a parametric study is carried out in this section. In order to find the individual impact of the design parameters, we kept all of the physical parameter values of the whole system constant (consistent with the values in Tables 1, 3, 4, and 5) except the parameter that we intend to evaluate its sole effect on the system. By using Equations (56) and (57), the expected value of the harvested electrical power may be obtained by the following equations.
Utilizing Equation (59), the mean squares of the normal stress components of the beams can be derived as It has to be clarified that since expected maximum stress is zero (i.e., E[σ max (t)] = 0), it cannot be utilized as a proper criterion to investigate the maximum deflection. Instead, since E σ 2 max (t) is of the order of magnitude of σ max (t), one can use it to provide an acceptable approximation of the order of magnitude of the maximum stress in the beams. So, the mean square stresses given in Equations (63) and (64) can be employed as a criterion to determine the physical failure limits during the random excitation of the system.

Parametric Study
In order to investigate the effects of different design parameters on the harvested electrical power, a parametric study is carried out in this section. In order to find the individual impact of the design parameters, we kept all of the physical parameter values of the whole system constant (consistent with the values in Tables 1 and 3-5) except the parameter that we intend to evaluate its sole effect on the system. In Figure 10, E P 1 (t) and E P 2 (t) have been plotted against the thickness ratio of the beams, s . In the simulations, h s . According to the figure, the expected values have major and minor peaks. As observed, the maximum of E P 1 (t) and E P 2 (t) occurs   Figure 10, which occurs between these two thickness ratios, is due to the phenomenon of mode veering. In this phenomenon, two modes of the dynamic system which are sufficiently close, excite and couple each other. This can lead to a higher vibration amplitude and greater harvested electrical power.
Similarly, the local maximum of the electrical power mean values, illustrated in Figure 10, which occurs between the thickness ratios of 0.2 and 0.4 is, again, due to the mode veering phenomenon when the first and second natural frequencies of the system become close at h  To justify the variation of the expected values presented above, the dependence of the first four natural frequencies on h (2) s /h (1) s is presented in Figure 11. As shown in the figure, the third and fourth natural frequencies approach closely to each other at the thickness ratio h  Figure 10, which occurs between these two thickness ratios, is due to the phenomenon of mode veering. In this phenomenon, two modes of the dynamic system which are sufficiently close, excite and couple each other. This can lead to a higher vibration amplitude and greater harvested electrical power.
Aerospace 2020, 7, x FOR PEER REVIEW 19 of 28 Figure 11. Normalized natural frequencies of the system versus the thickness ratio ℎ ( ) ℎ ( ) . Figure 12 shows the variation of ( ) and ( ) versus the stiffness of the interconnecting spring. It is easily observed that by increasing the rigidity of the beam coupling via increasing the spring stiffness, the mean value of the harvested power decreases. In other words, taking advantage of more compliant coupling between the beams, they undergo vibration with higher amplitude which leads to higher output power. Similarly, the local maximum of the electrical power mean values, illustrated in Figure 10, which occurs between the thickness ratios of 0.2 and 0.4 is, again, due to the mode veering phenomenon when the first and second natural frequencies of the system become close at h  Figure 12 shows the variation of E P 1 (t) and E P 2 (t) versus the stiffness of the interconnecting spring. It is easily observed that by increasing the rigidity of the beam coupling via increasing the spring stiffness, the mean value of the harvested power decreases. In other words, taking advantage of more compliant coupling between the beams, they undergo vibration with higher amplitude which leads to higher output power.  Figure 12 shows the variation of ( ) and ( ) versus the stiffness of the interconnecting spring. It is easily observed that by increasing the rigidity of the beam coupling via increasing the spring stiffness, the mean value of the harvested power decreases. In other words, taking advantage of more compliant coupling between the beams, they undergo vibration with higher amplitude which leads to higher output power. Another factor which plays a significant role in the efficiency of the coupled harvester, is the thickness of the piezoelectric layers. Figure 13 shows the dependence of the harvested power on this parameter. The mean power has an increasing-decreasing nature with respect to ℎ such that an optimum exists for each mean value. For instance, the maximum value of ( ) + ( ) for = 0.2 MΩ occurs at ℎ = 0.3348 mm. Another factor which plays a significant role in the efficiency of the coupled harvester, is the thickness of the piezoelectric layers. Figure 13 shows the dependence of the harvested power on this parameter. The mean power has an increasing-decreasing nature with respect to h p such that an optimum h p exists for each mean value. For instance, the maximum value of E P 1 (t) + E P 2 (t) for R l = 0.2 MΩ occurs at h p = 0.3348 mm. To determine why the maximum harvested power occurs at some intermediate value of h , we have to consider various phenomena. First, by decreasing h , the stress developed in the piezoelectric layer would decrease which leads to smaller electrical harvested power. However, as Equation (29) suggests, the generated voltages and, as a result, the harvested power depends on the To determine why the maximum harvested power occurs at some intermediate value of h p , we have to consider various phenomena. First, by decreasing h p , the stress developed in the piezoelectric layer would decrease which leads to smaller electrical harvested power. However, as Equation (29) suggests, the generated voltages and, as a result, the harvested power depends on the parameters ζ and ϑ k . Noting Equation (30), by increasing h p , ζ is increased which (since ζ acts similar to a damping factor in Equation (29)) is followed by a decrease in the harvested voltage. On the other hand, as shown in Figure 14, increasing ζ up to h p = 1.52 mm, leads to an increase in the coupling factor ϑ 2 , consequently leading to higher harvested voltages (see Figure 14). To determine why the maximum harvested power occurs at some intermediate value of h , we have to consider various phenomena. First, by decreasing h , the stress developed in the piezoelectric layer would decrease which leads to smaller electrical harvested power. However, as Equation (29) suggests, the generated voltages and, as a result, the harvested power depends on the parameters and . Noting Equation (30), by increasing h , is increased which (since acts similar to a damping factor in Equation (29)) is followed by a decrease in the harvested voltage. On the other hand, as shown in Figure (14), increasing ζ up to ℎ = 1.52 mm , leads to an increase in the coupling factor , consequently leading to higher harvested voltages.   Figure 15 demonstrates the expected value of the generated electrical power in both piezoelectric layers, for the case in which the electrical resistance, R l , has an optimal value of 0.5051 MΩ. It can be observed that by changing the electrical resistance from the mentioned nominal value to the current optimum value, the thickness of the piezoelectric layer h p , at which the mean of the total electric power is maximized, becomes 0.453 mm.
Aerospace 2020, 7, x FOR PEER REVIEW 21 of 28 Figure 15 demonstrates the expected value of the generated electrical power in both piezoelectric layers, for the case in which the electrical resistance, R , has an optimal value of 0.5051 MΩ . It can be observed that by changing the electrical resistance from the mentioned nominal value to the current optimum value, the thickness of the piezoelectric layer h , at which the mean of the total electric power is maximized, becomes 0.453 mm. Many other parameters can be used as design tools to increase the performance of the harvester. Among them are , and . In the next section, these parameters will be employed as design parameters in a genetic algorithm to find appropriate design values which maximize the total harvested power.

Design Optimization
The optimal selection of the geometrical and physical parameters of the energy harvester is critical in attaining high energy from the harvester. For example, as observed in the previous section, Many other parameters can be used as design tools to increase the performance of the harvester. Among them are R l , b s and b p . In the next section, these parameters will be employed as design parameters in a genetic algorithm to find appropriate design values which maximize the total harvested power.

Design Optimization
The optimal selection of the geometrical and physical parameters of the energy harvester is critical in attaining high energy from the harvester. For example, as observed in the previous section, a certain value of piezoelectric thickness can maximize the expected value of the harvested power. In the following study, the length of the beams and the piezoelectric patches are not considered as design parameters and their values are selected as those reported in Tables 1 and 3. Please note that in this table, L p < 0.1L s , which prevents charge cancellation when multiple structural modes are excited [27]. Moreover, to reduce the computational cost corresponding to the optimization procedure, the thickness of the thin beam is also considered to be the same as the one provided in Table 1. Assuming this to be the case, a physical understanding of the system reveals that under identical excitation conditions, the efficiency of the energy harvesting of the coupled system mainly depends on the thickness ratio of the beams, the stiffness of the interconnecting spring, the thickness of the piezoelectric layers, the width of the beams and the piezoelectric patches, as well as the load resistance. So, the following fitness function is considered.
In maximizing the fitness function, some constraints are taken into account to avoid an unrealistic physical system and violation of the modelling assumptions made throughout this paper. These constraints are summarized as: Imposing these constraints in Equation (66) prevents negative values for the thickness of either beam. It also guarantees that the upper beam remains thicker than the lower beam during the optimization process. The constraint given in Equation (67) can be restated as h (1) s ≤ 0.1L s and has been provided to preserve the Euler-Bernoulli beam assumptions. The upper limit in Equation (68) is devised to avoid a perfectly rigid connection between the two beams which in turn would lead to very high stiffness of the system. Moreover, the lower limit in this equation is provided to guarantee the coupling between the dynamics of the beams. The lower limits in Equations (69)-(71) have been selected to prevent negative values for the physical parameters. The upper limit in Equation (69) guarantees the satisfaction of the Euler-Bernoulli beam assumption for the piezo patches. In Equations (70) and (71), the width of the piezoelectric patches and that of the structures is forced to be less than 20% of their length, otherwise the beam assumption would be undermined, and a plate theory should be employed for system modelling. Furthermore, constraints in Equation (72) are imposed to provide the load circuit with arbitrary resistance from open to closed circuit. Finally, the constraints of Equations (73) and (74) have been considered to prevent yielding of the two beams while they are excited randomly. In the current problem, the yield stress of the both beams is assumed to be 250 MPa and a safety factor of two has been considered for the design of the system.
As shown earlier, the dynamic response of the system with characteristics given in Tables 1 and 3, can be accurately studied by considering only four modes. However, since the specifications of the system may alter during the optimization process, the number of effective modes may be affected. In order to make sure that sufficient modes are considered, a convergence study is carried out. To do so, in each iteration of the optimization process, the number of contributing modes n, is increased until the condition F n − F n+1 /F n ≤ 0.01 is satisfied. By employing the genetic algorithm (GA), the fitness function given in Equation (65) is maximized under constraints from Equations (66)-(73). The MATLAB default values, related to 'ga' function, are chosen for the genetic algorithm (GA). Simulation results, which are compiled in Table 8, demonstrate that in the optimum case, five modes participate in the dynamic response. It must be noted that the optimized characteristics of the system certainly depend on the fixed values chosen for the length of the piezoelectric, the length of the beam, and the thickness of the thin beam. If the harvester, with the parameters given in Tables 1 and 3 is replaced by the optimal harvester reported in Table 8, a remarkable improvement of 639.88% in the expected total harvested power is achieved. As is evident from Table 8, the optimal value related to the stiffness of the coupling spring is 665.94 N/m. At first glance, one may conclude that increasing the power of coupling by enhancing the spring stiffness from its lower band, 240 N/m, to the optimal value, 665.94 N/m, has increased the expected power obtained from the system, but this is not correct. If we compare the overall expected mean of the total power when k = 240 N/m with the expected power of the optimized system, we can see that the difference is negligible. This small difference is observed even when we put k = 0 instead of 240 N/m. This issue comes from the inherent properties of the GA since this process is a zero-order optimization algorithm. In the cases that the optimal values occur at the vertices of the region, the GA often requires a large number of iterations and normally convergence is assumed once the variation in optimal values of objective function falls below a threshold.
As discussed above, the coupling does not have a positive effect on enhancing the total harvested power of the optimized system which can be observed from Figure 12. This figure illustrates that when all of the other design variables are constant, increasing the coupling power decreases the mean expected harvested power of the system. Moreover, the optimum solution seems to have beams with very different thicknesses and therefore very little coupling in the modes. Consequently, there is not a significant advantage in coupling the two beams via a linear spring. This finding is opposite to what has been previously reported for energy harvesters under deterministic loads [12,13] since, in this state of excitation, changing the power of the coupling can tune one of the natural frequencies of the system to approach the single excitation frequency.
It is worth mentioning that the mode shapes of the optimized system do not have any node in the interval at which the piezo patches are attached. This verifies that charge cancellation does not happen in the system.

Conclusions
The importance of designing optimal energy harvesters to acquire green energy from environmental vibrations is well recognized. However, due to the random nature of ambient excitations, the electromechanical coupling of the governing equations, the multi-mode response of the system, and the interconnection of different bodies in the energy scavenging system, modelling and simulation of such harvesters remain complicated and non-trivial. A well-accepted approach to maximize the performance of the energy harvesting system is to match the resonant frequencies of the harvester with the frequency band of the excitation. To achieve this match, the design of energy harvesters should use different structural elements, mechanically coupled to each other to increase the number of natural frequencies in a pre-specified frequency band. In the current paper, the idea of designing an optimal structurally coupled energy harvester under random excitations was investigated. Coupled electromechanical equations of motion were derived using energy-based techniques. Then, using stochastic theories, the effect of random base motion on the harvested energy was studied analytically, and closed-form expressions were derived for the expected harvested power. Finally, a genetic algorithm was employed to determine the optimum physical and geometrical parameters which enhance the scavenged electrical power. The final optimized model demonstrated that the coupling of the structures lacks the potential to significantly improve the energy harvesting under random excitation. This finding is in stark contrast to previous studies on energy harvesting from deterministic excitations. The idea and the approach presented in this paper may be further employed to design other types of piezoelectric-based energy harvesters with complex structures undergoing random excitation. (1) and (2), eliminating exp( j * ω i t + φ i ) from both sides, and solving the resulting equations, one obtains ϕ (1) i (x) = C (1) i sin β i x + C (2) i cos β i x + C (3) i sinhβ i x + C (4) i cosh β i x (A2) ϕ (2) i (x) = C (5) i sin α i x + C (6) i cos α i x + C

By substituting Equation (A1), into Equations
i sinhα i x + C (8) i cosh α i x (A3) In these equations, C i , i ∈ N, j = 1, 2, . . . , 8, are constants that will be determined by satisfying boundary conditions. In addition, β i and α i are defined as β i = √ ω i and α i = 4 A Satisfying these boundary conditions leads to a set of linear algebraic homogenous equations of the form and A (i) is the coefficient matrix. By setting the determinant of the matrix A (i) to zero, the natural frequencies are obtained. Then the eigenvectors which characterize the mode shapes can be calculated.

Appendix B. Determination of the Location of the Neutral Surface
To determine the location of the neutral surface, the cross-sectional view of the beam and piezoelectric shown in Figure A1 is considered. In this figure, h sa is the distance between the lower surface of the substructure and the neutral axis, while h pa is the distance between the neutral axis and the upper surface of the piezoelectric.  Since there is no axial force, the following equation should be satisfied.
Using Hook's law, Equation (B-1) can be expressed in terms of strain as