Role of Electromechanical Coupling, Locomotion Type and Damping on the Effectiveness of Fish-Like Robot Energy Harvesters

Abstract

In this work, an investigation into the influence of prescribed motion on a body caudal fin aquatic unmanned vehicle (AUV) energy harvester is carried out. The undulatory–oscillation locomotion inspired by fishes actuates a composite beam representative of a spinal column with a piezoelectric patch. Two patch configurations—one at the head and tail—are considered for the AUV energy harvester, with a length that would not activate a harmonic in the system. An electromechanical model which accounts for the strain of the prescribed motion and the induced relative strain is developed. Discretizing the relative strain using Galerkin’s method requires a convergence study in which the impacts of the prescribed motion, dependent on the undulation and envelope of the motion, are investigated. The combination of prescribed motion and structural terms leads to a coupling that requires multiple investigations. The removal of the undulation of the system produces a more consistent response. The performances of the two different patch configurations undergoing different prescribed motions are studied in terms of coupled damping and frequency effects. An uncoupled Gauss law-based model is adopted to compare the performance of our approach and that of the coupled electromechanical model harvester. It is demonstrated that there is a complex interaction of the phases of the prescribed and relative motions of the structure which can lead to the development or destruction of the response of the total motion or voltage for the system. The results show that the structural damping and type of locomotion are the most influential parameters on the validity of the uncoupled approach. It is also found that the optimal resistances for the coupled and uncoupled representations are the same for the two motions and patch configurations considered.

1. Introduction

Aquatic unmanned vehicles (AUVs) are being developed and advancing into bioinspired body undulating–oscillatory locomotion, replicating fishes [1,2]. Typical aquatic systems are rigid-bodied, controlled with propellers and often tethered [3]. The advantage of bioinspiration is the replication of efficient swimming in close-quarter environments which are less disruptive than rigid propeller-driven systems [3]. The inclusion of soft and rigid materials helps to realize the structures necessary to generate the desired bioinspired locomotion [4,5]. The body shapes of these bioinspired AUVs are akin to the animals from which their motion is inspired [4,6]. These bioinspired body shapes are integrated with sensors that help monitor their surroundings [4,7,8]. The combination of body shape, actuation and sensors allows for an effective system [4].

Many species of animals have been used as sources of inspiration. To delineate this diverse range of animals in an orderly manner, a locomotion categorization has been proposed by Sfakiotakis et al. [9] and Salazar et al. [10]. A simplified version of this categorization is shown in Figure 1. Emphasis is placed on the fin oscillation caudal fin category because this is the most utilized locomotion type used in bioinspired AUVs [10]. Anguilliform locomotion represents the highly flexible-bodied AUVs that perform an eel-like swimming action. Subcarangiforms are less flexible-bodied than the Anguilliforms but allow for faster swimming predators. Carangiforms have a more restricted body motion and flexibility compared to Anguilliforms and Subcarangiforms, meaning that these animals are placed higher in the food chain. Thunniforms have very focused body motion to the peduncle, which creates a very powerful swimming motion for these top predators. The body composition and structure of individual species allow for the unique characteristics that cater to their survival.

More capable fish-like AUVs would have applications in civilian and military missions. Fish-like robots would be able to complete missions that require close-quarter maneuvering to monitor delicate or critical structures in coral reef research and oil field infrastructure [11]. Additionally, the more natural formation of these robotic systems would present the potential for them to integrate more naturally with their biological surroundings, offering valuable insights into undisrupted behavior [12]. The utility of systems that perform well in a challenging environment, as well as blending into the environment, offers a unique platform for a wide range of missions. The challenge of robotic systems operating in remote environments is the endurance necessary to complete the desired mission.

The endurance challenge of these systems requires a unique solution to increase the range of these aquatic systems. The flexibility and body actuation of bioinspired AUVs presents a potential platform for regenerative energy harvesting using an integrated energy harvesting patch. Works by Cha et al. [13] and Salazar et al. [14,15] have attempted to investigate this complex problem using a piezoelectric material. Piezoelectric transduction patches are becoming more commercially viable to be integrated into mechanical structures such as those found in AUVs. These patches are used to create an electrical potential when subjected to strain and have the capacity to power sensors that offset the power requirements of a primary battery. Cha et al. [13] constructed a physical system with a bimorph patch peduncle tail as a tapered cantilever beam, which was attached to a rotational spring to connect with the rotational excitation on the joint of this linkage. This induced bending in the patches on either side of the peduncle beam. Their study developed a coupled model for this structure and included effects from the aqueous environment. Salazar et al. [14,15] stated that the desired energy harvester should be excited by a defined spatio-temporal function that an AUV’s spinal column would move with if biomimetic in nature. In these works, only the spatio-temporal function was influenced by the electrical Gauss law equation. It was assumed that the structure would have very high damping and effectively eliminate the mechanical response to the prescribed motion; thus, the Gauss law equation would be dependent solely on the prescribed motion.

In this work, the development of an electromechanical piezoelectric energy harvester integrated into a representative AUV model is proposed in order to study the effect of damping, the type of locomotion and electromechanical coupling on the performance of the system. A comparative study between the coupled and uncoupled representations is carried out in order to determine the limits of applicability of the uncoupled Gauss law representation and to determine the interaction between the prescribed and relative motions of the fish-like robot energy harvesting system. The rest of this work is organized as follows. Section 2 demonstrates the representative model of the Carangiform motion that an AUV would be ideally using to swim. Section 3 breaks down the developed piezoelectric energy harvester model, showing the coupling between the relative and prescribed motions. Section 4 details the effect of the piezoelectric patch placement on the natural frequencies and mode shapes. Section 5 presents the results and discusses the phenomena from the perspective of how the interaction between the prescribed and relative motions of the AUV influence the effectiveness of the system for different kinds of undulatory–oscillatory representations. A particular focus is placed on how the structural damping affects the applicability of the uncoupled modeling. Section 6 concludes the work.

2. Bioinspired Prescribed Motion of Body Caudal Fin Carangiform AUV

AUVs that are developed with bioinspiration in mind commonly focus on the fin oscillation caudal fin category, as fishes in this category are proficient and very capable swimmers. The bodies of these animals participate in the kind of locomotion known as body caudal fin (BCF) swimming. These animals have a skeleton that supports a musculature that controls excitation from the head to the tail of the animal along a spinal chain. The spinal chain allows for flexible yet strong support for the body to perform an undulating–oscillatory motion. The motion of the spinal chain is dependent on the species and the architecture of these animals’ internal structures. Indeed, the variability of the locomotion categorization is related to the number of segmentations of the spinal chain, where Anguilliforms tend to have many and Thunniforms have far fewer segmentation [4,10]. These animals engage their musculature and translate a moving wave along the body to the tail tip to produce thrust for swimming. BCF locomotion can be represented as a sinusoidal waveform which is restricted by a polynomial envelope, as shown in Equation (1). The respective constraints of this motion are dependent on the animal to be replicated [15].

$$v_p\left(x,t\right)=L\left(a_0+a_1\left(\frac{x}{L}\right)+a_2{\left(\frac{x}{L}\right)}^2\right)\cos \left(\frac{Kx}{L}-{\omega }_{p}t\right)$$

(1)

The spatio-temporal function, $v_p\left(x,t\right)$, is the prescribed motion which a bioinspired AUV would be trying to replicate. The constants $a_0$, $a_1$ and $a_2$ determine the constraining envelope of the prescribed motion, given as 0.02, −0.0825, and 0.1625 for Carangiform motion, respectively [16]. The flexibility of the spinal column replicating a periodic wave is controlled by the undulatory term K, where one periodic wave equates to this undulating term being $2\pi /L$. This motion has been shown to be capable in AUV designs, and many systems are being developed to perform this swimming action. AUVs can have rigid segmented chains or flexible soft structures. The sinusoidal function is manipulated in accordance with the constraint that only one period should occur over the length, as the spinal column cannot be too flexible if it is to remain within these constraints of AUV functional design. The frequency ($f_p$) of oscillation is defined to be constrained between 2–5 Hz; thus ${\omega }_p=$ 12.566–62.832 rad/s. Figure 2 demonstrates multiple time instants of a defined Carangiform motion with the described envelope terms and one period of spinal undulation.

3. Piezoelectric Energy Harvester Methods

The development of piezoelectric energy harvesters has reached a transition point, leading to possible applications [17,18,19]. Lead zirconium-titanate (PZT) is considered to be one of the best-performing piezoelectric materials and is suitable for a wide range of applications [17]. PZT materials depend on the organization of precursor elements into the desired perovskite structure [20]. The perovskite structure can be doped with additional elements to boost the percentage of this structure or influence the properties of the material for better performance [21]. The perovskite structure is organized into larger bulk transduction materials that are integrated into patches of two different structures: thin-film or macro-fiber composite (MFC). These different patch architectures give a variability of patch flexibility and performance that is suitable for unique application scenarios depending on the excitation. The MFC developed by Smart Material Corporation is highly flexible due to its composite nature and is an ideal selection for excitations that subject these patches to large strain. Research works that define the properties of materials in piezoelectric energy harvesters can be found in [22,23,24,25,26,27,28]. Next, the coupled and uncoupled representations for the determination of the levels of the harvested power from the prescribed motion of fish-like robots are presented.

3.1. Electromechanical Coupled Piezoelectric AUV Energy Harvester Modeling

The application in a bioinspired AUV energy harvester of a piezoelectric patch on a spinal column that is flexible has been investigated in [14]. It was demonstrated in [13] that the optimal piezoelectric patch length is independent of the scale of the AUV. Additionally, for Carangiform motion, the optimal patch placement is at the location of the greatest strain: the tail tip. The authors investigated an AUV energy harvester under the assumption that the shunt damping is negligible due to the huge structural damping in the system’s design; thus, purely the energy generation is described. However, when neglecting the natural relative motion of the structure, possible unexpected and damaging phenomena of forced motion could be overlooked. Therefore, this work investigates the coupling of a piezoelectric transduction patch application on an AUV backbone substrate subjected to multiple undulating–oscillating motions to determine the possible occurrences in a bioinspired AUV. The structure of the proposed AUV energy harvester is shown in Figure 3.

This body caudal fin energy harvester (BCFEH) attached at the tail tip undergoes the prescribed Carangiform motion ($v_p$), as shown in Figure 3b. The natural frequencies of the BCFEH are dependent on the materials and the structural configuration. Many of the fixed dimensions are shown in

Table 1. Fixed characteristics of the body caudal fin energy harvester (BCFEH) adopted from [15].

Symbol	Parameter	Value
$b_p$	Width of piezoelectric layer (mm)	7
$b_s$	Width of substrate layer (mm)	7
$h_p$	Thickness of piezoelectric layer (mm)	0.2
$h_s$	Thickness of aluminum substrate (mm)	0.6
$E_p$	Piezoelectric material Young’s modulus (GN/m2)	30.336
$E_s$	Aluminum Young’s modulus (GN/m2)	69.5
${\rho }_s$	Density of aluminum substrate (kg/cm3)	2700
${\rho }_p$	Density of piezoelectric patch (kg/cm3)	5440
$m_1$	$\mathrm{Mass} \mathrm{per} \mathrm{unit} \mathrm{length} (0\le x\le L_1$)	${\rho }_s{b}_s{h}_s$
$m_2$	$\mathrm{Mass} \mathrm{per} \mathrm{unit} \mathrm{length} (L_1<x\le L_2$)	${\rho }_s{b}_s{h}_s+{\rho }_p{b}_p{h}_p$
$\overline{y}$	Position to the neutral axis	$\frac{E_p{h}_p\left(h_p+h_s\right)}{2\left(E_p{h}_p+E_s{h}_s\right)}+\frac{h_s}{2}$
$y_1$	Position of aluminum relative to the neutral axis	$h_s-\overline{y}$
$y_2$	Position of piezoelectric layer relative to the neutral axis	$(h_s+h_p)-\overline{y}$
$y_0$	Start of aluminum layer relative to neutral axis	$-\overline{y}$
$e_{31}$	$\mathrm{Piezoelectric} \mathrm{stress} \mathrm{constant} (\mathrm{C}/{\mathrm{m}}^2$)	$-5.16$
${ϵ}_{33}^{\sigma }$	$\mathrm{Permittivity} \mathrm{constant} (\mathrm{nF}/\mathrm{m}$)	12.563

. Figure 3 shows the patch at the tail tip, in accordance with the previous work; however, a patch at the head is also considered to compare how the electromechanical coupling could influence the energy harvester’s performance. The length of the patch is defined as 0.1 m for both cases to avoid the resonance region and resulting high-amplitude oscillations. The aim is to consider an off-resonance scenario and study the impacts of the undulatory–oscillatory motion and electromechanical coupling on the effectiveness of the system and the validity of the uncoupled modeling when low damping takes place. The backbone is considered as a cantilever beam with a clamped end condition at $L_0=0 \mathrm{m}$.

The response of the piezoelectric materials undergoing stress has commonly been accounted for in the Gibb’s free energy (G), such as in works by Amin et al. [29], Sherrit and Mukherjee [30] and Wang et al. [31]. This formulation can be expanded to higher-order piezoelectric, permittivity and electrostrictive terms, such as in the works of Guyomar et al. [32,33] and Joshi [34], to be

$$G=-e_{31}{\epsilon }_{11}{E}_3-\frac{1}{2}{\in }_{33}^{{\sigma }_c}{E}_3^2+\frac{1}{2}\beta {\epsilon }_{11}{E}_3^2-\frac{1}{2}\gamma {\epsilon }_{11}^2{E}_3-\frac{1}{6}\delta E_3^3$$

(2)

where G depends on the piezoelectric strain constant ($e_{31}$), the strain of the material is subjected to (${\epsilon }_{11}$), the electric field ($E_3$), permittivity under constant stress (${\in }_{33}^{{\sigma }_c}$), nonlinear electrostriction constant ($\beta$), nonlinear piezoelectricity constant ($\gamma$) and nonlinear permittivity constant ($\delta$). In this study, to focus on the electromechanical coupling effects on the harvester’s response and its comparison to uncoupled modeling, the higher-order nonlinear terms are neglected, and thus only conventional linear piezoelectric and permittivity terms are considered as follows:

$$G=-e_{31}{\epsilon }_{11}{E}_3-\frac{1}{2}{\in }_{33}^{{\sigma }_c}{E}_3^2$$

(3)

The electric displacement ($D_3$) is found by deriving the function of Equation (3) with respect to the electric displacement [34]:

$$D_3=-\frac{\partial G}{\partial E_3}$$

(4)

which yields [35]

$$D_3=e_{31}{\epsilon }_{11}+{\in }_{33}^{{\sigma }_c}{E}_3$$

(5)

The stress of the piezoelectric is determined by deriving Equation (3) with respect to strain, as shown in Equation (6).

$${\sigma }_{11}^p=E_p{\epsilon }_{11}+\frac{\partial G}{\partial {\epsilon }_{11}}$$

(6)

Thus, the piezoelectric material’s stress response is expressed as [36]

$${\sigma }_{11}^p=E_p{\epsilon }_{11}-e_{31}{E}_3$$

(7)

The strain of the BCFEH backbone is determined to be dependent on the conventional Euler–Bernoulli beam theory, as the backbone structure is comparable to a long, thin beam. Forcing the prescribed fish-like motion onto the structure and incorporating the response into the strain, Equation (8) is determined as

$${\epsilon }_{11}=-y\left(v_r^{″}+v_p^{″}\right)$$

(8)

where y is the relationship of the backbone and patch thickness to the structural neutral axis.

To accurately account for the multiple excitations and response of the structure, an energy model formulation is adopted from [15,37,38]. The kinetic energy (T), shown in Equation (9), of the BCFEH is dependent on ${\dot{v}}_r$ and ${\dot{v}}_p$, which are the velocity of the relative and prescribed motion, respectively. Equation (10) is the potential energy ($\mathsf{\Pi }$) of the BCFEH as it undergoes motion.

$$T=\frac{1}{2}m\int {\left(\dot{v_r}+{\dot{v}}_p\right)}^{2}dx$$

(9)

$$\mathsf{\Pi }=\frac{1}{2}\left[\iint {\sigma }_{11}^s{\epsilon }_{11}d{A}_{s}dx+\iint {\sigma }_{11}^p{\epsilon }_{11}d{A}_{p}dx-\iint E_3{D}_{3}d{A}_{p}dx\right]$$

(10)

The mass, m, is dependent on the location along the length, where $m=m_1$ for the first section of the backbone without the piezoelectric patch and $m=m_2$ for the tail section covered by the piezoelectric patch. The stress of the structure (${\sigma }_{11}^s$) is related through Hooke’s law. The electric field is related to the voltage generation through the piezoelectric patch thickness by

$$E_3=-\frac{V\left(t\right)}{h_p}$$

(11)

The relative displacement from the response of the structure is accounted for in a Galerkin discretization as

$$v_r\left(x,t\right)={\displaystyle \sum }_{i=1}^N{\varphi }_{ji}\left(x\right)q_i\left(t\right)$$

(12)

The displacement is dependent on the mode shape (${\varphi }_{ji}\left(x\right)$), where subscripts are dependent on both the section of the beam and mode of vibration; $j=1$ or $2$ corresponds to the first or second section of the beam’s length, respectively, and $i$ denotes the vibration mode number. The number of modes used is determined by N. Modal coordinates are denoted as $q_i\left(t\right)$. The natural frequencies and mode shapes are determined for a cantilever beam with different layer regions by an eigenvalue system of equations considering typical boundary conditions that satisfy the orthogonality conditions for the systems’ natural frequencies [22].

Equations (5), (7) and (8) are combined into Equations (9) and (10) with the consideration of Equations (11) and (12). The electromechanical equations of motion are found using the Euler–Lagrange principle as follows:

$$\frac{d}{dt}\left(\frac{\partial L}{\partial {\dot{q}}_i}\right)-\frac{\partial L}{\partial q_i}=\frac{\delta W_{nc}}{\delta q_i}=-2\xi {\omega }_i{\dot{q}}_i$$

(13)

$$\frac{d}{dt}\left(\frac{\partial L}{\partial \dot{\lambda }}\right)-\frac{\partial L}{\partial \lambda }=\frac{\delta W_{nc}}{\delta \lambda }=\frac{-V}{R}$$

(14)

where ($L=T-\mathsf{\Pi }$) and $V\left(t\right)=\dot{\lambda }$. The damping ratio ($\xi$) is defined for the structure and R denotes the load resistance of the circuit.

Using the Euler-Lagrange equations, the following reduced-order model is obtained:

$${\ddot{q}}_i+2\xi {\omega }_i\dot{q_i}+{\omega }_i^2{q}_i-\theta \left\{{\varphi }_i^{\prime }\left(L_{x_f}\right)-{{\varphi }_i}^{\prime }\left(L_{x_i}\right)\right\}V\left(t\right)=-m_1{{\displaystyle \int }}_0^{L_1}\left({\ddot{v}}_p{\varphi }_{1i}\right)dx$$

$$-m_2{{\displaystyle \int }}_{L_1}^{L_2}\left({\ddot{v}}_p{\varphi }_{2i}\right)dx-E{I}_1{{\displaystyle \int }}_0^{L_1}\left\{v_p^{″}{\varphi }_{1i}^{″}\right\}dx-E{I}_2{{\displaystyle \int }}_{L_1}^{L_2}\left\{v_p^{″}{\varphi }_{2i}^{″}\right\}dx$$

(15)

$$C_p\dot{V}\left(t\right)+\frac{V\left(t\right)}{R}+\theta {{\displaystyle \int }}_{L_{x_i}}^{L_{x_f}}\left\{{\displaystyle \sum }_{i=1}^N{\varphi }_i^{″}\dot{q_{𝚤}}\right\}dx=-\theta {{\displaystyle \int }}_{L_{x_i}}^{L_{x_f}}\left\{\dot{v_p^{″}}\right\}dx$$

(16)

The piezoelectric coupling term and capacitance of the piezoelectric patch are defined as

$$\theta =\frac{e_{31}{b}_p\left(y_2+y_1\right)}{2}$$

$$C_p=\frac{{\in }_{33}^{{\sigma }_{ }}{b}_p\left(L_{x_f}-L_{x_i}\right)}{h_p}$$

It should be noted that Equation (15) is dependent on the mode shape and natural frequency number; thus, a system of ($N+1$) coupled equations is constructed. Moreover, it is clear that the prescribed motion plays an important role in the structurally dominant Equation (15). Indeed, the prescribed motion is a source of excitation and may result in the presence of high amplitudes of relative motion, and thus the level of generated voltage is affected by both relative and prescribed motions, as shown in Equation (16). It follows from Equation (15) that the prescribed motion results in two kinds of excitations for the fish-like robot from a structural point of view. The first excitation is due to the acceleration of the prescribed motion, and the second source of excitation is due to the undulatory strain effects of the prescribed motion. These two kinds of excitations are evaluated in the next sections, as well as how the type of undulatory–oscillatory motion of the fish-like robot system and structural damping affect the effectiveness of the energy harvesting system.

3.2. Uncoupled Energy Harvesting Modeling

Considering the high damping in the structural design of these types of fish-like robot is important in order to avoid any source of external excitations. Following the works presented in [13,14], an uncoupled modeling can be performed based on the Gauss law equation by assuming that the prescribed motion is the only source of excitation and that relative motion will not take place, and thus the electromechanical coupling is neglected. This means that Equation (15) is not needed, and the relative motion effects in Equation (16) do not exist. In other words, a comparative example of electrical potential harvested solely from the prescribed motion requires an assumption that there is no coupling from the structure. The possibility of having a structure with an internal mechanism with high damping could allow for such a case [39,40]. This would require a damping mechanism that could silence the induced natural mechanical motion from the prescribed motion but still allow the structure to move with the prescribed motion. As $\xi \to \infty$, the mechanical contribution in Equation (16) decreases until it becomes negligible, and the equation follows the uncoupled scenario with only the Gauss law equation [13].

The prescribed motion is the only excitation affecting the piezoelectric layer. Following the work presented in [13], the trigonometric terms are separated to extract the constant parameters $K_1$ and $K_2$.

$$\begin{gathered} K_1{}_{ }={\left.\left(-\left[\left(\frac{a_1}{L}+\frac{2a_2}{L}{\left(\frac{x}{L}\right)}^{ }\right)\right]\cos \left(Kx\right) \\ +\left[K\left(a_0+a_1\left(\frac{x}{L}\right)+a_2{\left(\frac{x}{L}\right)}^2\right)\right]\sin \left(Kx\right)\right)\right|}_{L_{x_i}}^{L_{x_f}} \end{gathered}$$

(17)

$$K_2{}_{ }={\left.\left(\left[\left(\frac{a_1}{L}+\frac{2a_2}{L}{\left(\frac{x}{L}\right)}^{ }\right)\right]\sin \left(Kx\right)+\left[K\left(a_0+a_1\left(\frac{x}{L}\right)+a_2{\left(\frac{x}{L}\right)}^2\right)\right]\cos \left(Kx\right)\right)\right|}_{L_{x_i}}^{L_{x_f}}$$

(18)

A closed-form solution is employed to obtain the steady-state solution. Thus, the final expression of the generated voltage for the aquatic unmanned fish energy harvester is given by

$$V\left(t\right)=\left|\theta \right|{\omega }_{p}R\sqrt{\frac{K_1^2+K_2^2}{1+{\left(C_p{\omega }_{p}R\right)}^2}}\sin \left({\omega }_{p}t+\varphi \right)$$

(19)

The root mean squared for the generated voltage is obtained for a sinusoidal as

$$V_{RMS}=\frac{\left|\theta \right|{\omega }_{p}R}{\sqrt{2}}\sqrt{\frac{K_1^2+K_2^2}{1+{\left(C_p{\omega }_{p}R\right)}^2}}$$

(20)

The average power is then determined as

$$P_{avg}=\frac{{\left[V_{RMS}\right]}^2}{R}=\frac{{\theta }^2{\omega }_p^{2}R\left(K_1^2+K_2^2\right)}{2(1+{\left(C_p{\omega }_{p}R\right)}^2}$$

(21)

The uncoupled $V_{RMS}$ and $P_{avg}$ in Equations (20) and (21) are compared to the coupled system found through the Euler–Lagrange methodology.

4. Effect of Patch Placement on Linear Characteristics of the System

The development of the BCFEH requires an understanding of how the various excitations can affect the response of the fish-like energy harvester. If the excitation of the prescribed motion or environmental response equals one of the natural frequencies of the harvester, structural damage could result. The natural frequencies and mode shapes of the BCFEH are determined based on the cantilever configuration with partial patch coverage at the head or the tail tip. The partial patch head configuration is considered to correspond to $L_{x_i}=0 \mathrm{m}$ and $L_{x_f}=0.03 \mathrm{m}$, and the locations for the tail configuration are $L_{x_i}=0.07 \mathrm{m}$ and $L_{xf}=0.1 \mathrm{m}$. The BCFEH is defined to operate with prescribed frequencies less than 10 Hz. To avoid the natural frequencies of the system causing resonance that could damage the system during the defined operating conditions, the BCFEH was defined to be 0.1 m in length, and a patch of 0.03 m in length was attached at the head and tail.

Table 2. First three natural frequencies for the different configurations.

L (m)	First Natural Frequency (Hz)	Second Natural Frequency (Hz)	Third Natural Frequency (Hz)
0.1	Head Attachment
	59.4099	338.3119	898.0767
0.1	Tail Attachment
	42.1036	303.1340	874.8842

lists the first three natural frequencies for the BCFEH with these two respective patch attachment locations to show that the energy harvester worked in a region that was not susceptible to resonance. Multiple modes were needed to perform the convergence analysis. Figure 4 shows the first seven mode shapes for the two patch configurations.

Figure 4a demonstrates that the head patch configuration decreased the non-normalized amplitude on the first section of the length when compared to the tail configuration shown in Figure 4b. This is due to the increased stiffness from the patch coverage at the head section. Alternatively, the head patch configuration had a larger non-normalized amplitude at the tail tip, as this section only included substrate and was therefore less stiff. The tail patch configuration exhibited the opposite phenomenon along the length as the non-normalized amplitude was larger near the head and reduced near the tail.

5. Importance of Prescribed Motion on Structural Motion and Limits to the Applicability of the Uncoupled Modeling: Role of Forced Actuation and Damping

The forced actuation that is desired for the BCFEH to perform missions is important to investigate as it creates more forcing terms in Equations (15) and (16). These additional terms depend on how the $v_p$ is defined.

Table 3. Classifications of $v_p$ cases.

Case	Equation	Description
A	$L\left(a_o+a_1\left(\frac{x}{L}\right)+a_2{\left(\frac{x}{L}\right)}^2\right)\cos \left(Kx-{\omega }_{p}t\right)$	Undulating–oscillatory motion with original backbone polynomial.
B	$L\left(a_o+a_1\left(\frac{x}{L}\right)\right)\cos \left(Kx-{\omega }_{p}t\right)$	Undulating–oscillatory motion with linear backbone constraint as $a_2=0$.
C	$L\left(a_o\right)\cos \left(Kx-{\omega }_{p}t\right)$	Undulating–oscillatory motion with constant backbone constraint as $a_2=a_1=0$.
D	$L\left(a_o+a_1\left(\frac{x}{L}\right)+a_2{\left(\frac{x}{L}\right)}^2\right)\cos \left({\omega }_{p}t\right)$	Oscillatory motion with original backbone polynomial.
E	$L\left(a_o+a_1\left(\frac{x}{L}\right)\right)\cos \left({\omega }_{p}t\right)$	Oscillatory motion with linear backbone envelope.
F	$L\left(a_o\right)\cos \left({\omega }_{p}t\right)$	Oscillatory motion with constant backbone constraint.

describes the different prescribed motion classifications that are considered: cases A–C are undulating–oscillatory motion in which the envelope of the motion is altered, and cases D–F are oscillatory motion in which the undulatory term K equals 0, and the envelope of the motion alters incrementally as in Cases A–C.

Figure 5 presents visualization of the cases of motion defined in Table 3 for multiple time instants. Case A is shown in Figure 5a, which is the undulation–oscillatory motion defined for Carangiform fishes. Envelope terms $a_0$, $a_1$ and $a_2$ and the undulatory term $K$ are defined based on this biological motion. Case B shown in Figure 5b occurs when $a_2=0$; thus, a linear envelope is achieved. Case C in Figure 5c is when $a_1=a_2=0$, creating a constant envelope constraint for the undulation–oscillatory motion. Figure 5d exhibits the oscillatory motion of Case D as $K=0$ and is constrained by the original envelope terms. Since this case has a quadratic envelope, the beam is still forced with bending but does not have the undulatory motion present in cases A–C. Figure 5e shows that case E has a linear envelope and that there is no bending actuation as in case D. Figure 5f shows case F and has a constant envelope within which the beam oscillates. Cases E and F would generate the more common base excitation forcing terms in Equation (15).

5.1. Convergence Analysis Investigation for Undulatory–Oscillatory Prescribed Motion

When discretizing the equations of motion, it is essential to investigate the convergence of the BCFEH response for both head and tail patch configurations while undergoing the prescribed motion cases in Table 3. The prescribed motion interacts with the relative structural motion and requires a defined number of modes to obtain a consistent response for the fish-like energy harvester. The maximum relative displacement, average harvested power and generated voltage responses are investigated for the coupled BCFEH considering $N=1\to 7$. Case A motion is considered first in the convergence analysis. Figure 6 and Figure 7 show the induced relative motion and harvested power for the head and tail patch configurations, respectively. For the excitation frequency, $f_p$, the range considered is 0.25–10 Hz, $\xi =0.3$ and $R={10}^5 \mathsf{\Omega }$. The time histories consider $f_p$ = 5 Hz and a time length of 15 s. The uncoupled average harvested power is also shown for comparison. It should be mentioned that only off-resonance scenarios are considered in order to check the validity of the uncoupled modeling versus low values of damping for the coupled system.

Figure 6a shows that the maximum $v_r\left(L\right)$ in response to the prescribed motion is two times larger than that of the defined $v_p\left(L\right)$ amplitude of 0.01 m. The case A patch attached at the head $v_r\left(L\right)$ does not converge to a single response. The $P_{avg}$ performance does reach more of a convergence with a higher number of modes, as depicted in Figure 6b. The time history shows that the number of modes affects not only the amplitude but also the phase of the response, as shown in Figure 6c,d. The coupled system exhibits higher amounts of harvested power than the uncoupled case for all numbers of modes.

Figure 7 shows the tail attachment configuration BCFEH response for case A. Regarding the head patch configuration, the relative displacement for this tail configuration under the prescribed motion in case A exhibits a larger amplitude than the input 0.01 m tail amplitude. The relative displacement and average power do not reach complete convergence. The $P_{avg}$ approaches a convergence at lower values than the uncoupled case. The time histories for the relative displacement and average harvested power show that the amplitude and phase approach convergence with a higher number of modes, as the head patch configuration did. However, the tail patch coupled system is more affected by the number of modes. Clearly, case A undulatory–oscillatory prescribed motion affects both forcing terms in Equation (15), which requires the reduced-order model to have a higher number of modes in the Galerkin discretization, although the system of equations is linear. It is clear that the activation of the second derivative of the prescribed motion in case A means that the forcing is highly dependent on the mode number. More investigation on this is shown in the following sections.

The convergence of cases B and C was also investigated in the same manner; however, the results were unrealistic and are omitted from this investigation. Amplitudes for the $v_r\left(L\right)$ and $P_{avg}$ were too large and were unrealistic. The combination of the mode shapes with the altered polynomial envelope produced terms that influenced the system too much. This unrealistic behavior is discussed below. Therefore, case A serves as an appropriate case for further investigation to determine the impacts of the $K$ undulatory term when $N=7$. Next, the investigation of case D is considered in Figure 8 and Figure 9 to determine the effects of prescribed motion on the system’s response when the undulatory term is considered to be zero. The values of $f_p$ and $\xi$ are considered to be the same as in case A. Time histories consider an $f_p$ value of 5 Hz and a time length required for a response to steady state which is found to be longer than for case D.

The convergence head patch configuration considering case D motion is shown in Figure 8. A very different response of the system to case A is evident. For this motion case, the relative displacement is much less affected and deviates slightly from a constant response for all prescribed frequencies considered. The average harvested power approaches a converged response that is higher than the uncoupled case. The phase of the $v_r\left(L\right)$ and voltage time history are not affected in the same manner as in case A. This demonstrates that the undulatory term, $K$, has a more complex influence on the system’s response. The tail patch configuration has a very similar response to the head patch case, as shown in Figure 9.

Cases E and F are representative of a typical base excitation. It is found that when the complexity of the undulation and the forced curvature of the polynomial envelope are removed, the system generates a converged response in a few modes due to the linear characteristics of the reduced-order model. Since the uncoupled power model requires at least a quadratic envelope constraint, the estimated response of the system is zero. Therefore, to compare the uncoupled system to a coupled case, case D is selected for further investigation.

Table 4. Table of convergence for case D: maximum $v_r\left(L\right)$ and $P_{avg}$, $f_p=5 \mathrm{Hz}$ and $\xi =0.3$.

	Head		Tail
	${\mathit{v}}_{\mathit{r}}\left(\mathit{L}\right)\times {10}^{-2} \left(\mathbf{m}\right)$	$\frac{\left|\mathit{D}\mathit{i}\mathit{f}\mathit{f}\right|}{\mathit{V}\mathit{a}{\mathit{l}}_{\mathit{N}=6}}$	${\mathit{v}}_{\mathit{r}}\left(\mathit{L}\right)\times {10}^{-2} \left(\mathbf{m}\right)$	$\frac{\left|\mathit{D}\mathit{i}\mathit{f}\mathit{f}\right|}{\mathit{V}\mathit{a}{\mathit{l}}_{\mathit{N}=6}}$
Case A $R={10}^5 \mathsf{\Omega }$
N = 6	2.1668	0.99%	2.2029	0.93%
N = 7	2.1453		2.1824
Case D $R={10}^5 \mathsf{\Omega }$
N = 6	1.6198	0.12%	1.6207	0.093%
N = 7	1.6217		1.6222
	Head		Tail
	$P_{avg}$ (W)	$\frac{\left|Diff\right|}{Va{l}_{N=6}}$	$P_{avg}$ (W)	$\frac{\left|Diff\right|}{Va{l}_{N=6}}$
Case A $R={10}^5 \mathsf{\Omega }$
N = 6	$3.3166\times {10}^{-3}$	2.62%	$4.5115\times {10}^{-3}$	11.45%
N = 7	$3.4035\times {10}^{-3}$		$3.9949\times {10}^{-3}$
Case D $R={10}^5 \mathsf{\Omega }$
N = 6	$3.3612\times {10}^{-4}$	3.68%	$3.4205\times {10}^{-4}$	1.49%
N = 7	$3.4849\times {10}^{-4}$		$3.4715\times {10}^{-4}$

gives values for $N$ of 6 and 7 in cases A and D for the two patch configurations. A single prescribed frequency of 5 Hz is selected. This table shows that the variation between these higher modes in case A is larger for $v_r\left(L\right)$ and $P_{avg}$ than case D. The $v_r\left(L\right)$ difference between the head and tail configurations is perceived to be because of the mode shape. The head configuration is more flexible in the tail region and could lead to a slightly larger difference in the relative displacement. The tail configuration is stiffer at the tail, and thus a more rigid system is created that converges with a smaller number of modes at this point. There is a definite impact of the undulation term K that can be seen in the case A $P_{avg}$ values. The case A tail patch configuration shows the greatest variation of harvested power out of all the cases investigated at 11.45%. When $K=0$, the impact on the harvested power variation for the higher modes is reduced. Further investigation is required to explain this phenomenon.

To examine the factors influencing the BCFEH, the forcing terms in Equation (15) are clarified in terms of their magnitude depending on the considered motion case and are detailed in

Table 5. Table of forcing terms depending on mode number for the head patch.

	$D_0=D_1\cos \left(\omega t\right)+D_{2}sin\left(\omega t\right)=Dcos\left(\omega t+{\varphi }_1\right)$; where $D=\sqrt{D_1^2+D_2^2}$ ${\mathsf{\Lambda }}_0={\mathsf{\Lambda }}_1{\omega }_p^{2}L\cos \left(\omega t\right)+{\mathsf{\Lambda }}_2{\omega }_p^{2}Lsin\left(\omega t\right)=\mathsf{\Lambda }{\omega }_p^{2}Lcos\left(\omega t+{\varphi }_2\right)$; where $\mathsf{\Lambda }=\sqrt{{\mathsf{\Lambda }}_1^2+{\mathsf{\Lambda }}_2^2}$
	Head
	Case A	Case B	Case C	Case D	Case E	Case F
$i=1$	$D=51.0491$ $\mathsf{\Lambda }=1.5124\times {10}^{-5}$	$D=27.6550$ $\mathsf{\Lambda }=2.5322\times {10}^{-5}$	$D=38.3298$ $\mathsf{\Lambda }=1.7342\times {10}^{-5}$	$D=34.9648$ $\mathsf{\Lambda }=1.3927\times {10}^{-4}$	$D=0$ $\mathsf{\Lambda }=1.0514\times {10}^{-4}$	$D=0$ $\mathsf{\Lambda }=5.2968\times {10}^{-5}$
$i=2$	$D=376.2971$ $\mathsf{\Lambda }=6.0085\times {10}^{-5}$	$D=349.5832$ $\mathsf{\Lambda }=5.9494\times {10}^{-5}$	$D=337.7887$ $\mathsf{\Lambda }=1.8576\times {10}^{-5}$	$D=66.2740$ $\mathsf{\Lambda }=1.2492\times {10}^{-5}$	$\mathsf{\Lambda }=2.4546\times {10}^{-6}$	$\mathsf{\Lambda }=3.4066\times {10}^{-5}$
$i=3$	$D=1.3211\times {10}^3$ $\mathsf{\Lambda }=8.6913\times {10}^{-5}$	$D=1.1278\times {10}^3$ $\mathsf{\Lambda }=5.4179\times {10}^{-5}$	$D=829.2084$ $\mathsf{\Lambda }=5.6216\times {10}^{-5}$	$D=143.2050$ $\mathsf{\Lambda }=1.4790\times {10}^{-5}$	$\mathsf{\Lambda }=1.0209\times {10}^{-5}$	$\mathsf{\Lambda }=2.1958\times {10}^{-5}$
$i=4$	$D=3.5329\times {10}^3$ $\mathsf{\Lambda }=1.7550\times {10}^{-5}$	$D=2.7855\times {10}^3$ $\mathsf{\Lambda }=2.7622\times {10}^{-5}$	$D=678.1077$ $\mathsf{\Lambda }=2.9229\times {10}^{-5}$	$D=256.9345$ $\mathsf{\Lambda }=7.6621\times {10}^{-6}$	$\mathsf{\Lambda }=9.8027\times {10}^{-6}$	$\mathsf{\Lambda }=1.5544\times {10}^{-5}$
$i=5$	$D=4.1794\times {10}^3$ $\mathsf{\Lambda }=2.1318\times {10}^{-5}$	$D=1.7432\times {10}^3$ $\mathsf{\Lambda }=1.3214\times {10}^{-5}$	$D=777.3526$ $\mathsf{\Lambda }=1.7981\times {10}^{-5}$	$D=167.5726$ $\mathsf{\Lambda }=9.3540\times {10}^{-6}$	$\mathsf{\Lambda }=8.8458\times {10}^{-6}$	$\mathsf{\Lambda }=1.2392\times {10}^{-5}$
$i=6$	$D=4.2146\times {10}^3$ $\mathsf{\Lambda }=7.8878\times {10}^{-6}$	$D=3.1910\times {10}^3$ $\mathsf{\Lambda }=1.941\times {10}^{-5}$	$D=967.2774$ $\mathsf{\Lambda }=1.5006\times {10}^{-5}$	$D=352.5435$ $\mathsf{\Lambda }=7.3019\times {10}^{-6}$	$\mathsf{\Lambda }=7.7865\times {10}^{-6}$	$\mathsf{\Lambda }=1.0187\times {10}^{-5}$
$i=7$	$D=5.5492\times {10}^3$ $\mathsf{\Lambda }=1.0522\times {10}^{-5}$	$D=2.9218\times {10}^3$ $\mathsf{\Lambda }=9.7966\times {10}^{-6}$	$D=721.7400$ $\mathsf{\Lambda }=1.0481\times {10}^{-5}$	$D=446.7825$ $\mathsf{\Lambda }=6.9745\times {10}^{-6}$	$\mathsf{\Lambda }=6.6603\times {10}^{-6}$	$\mathsf{\Lambda }=8.3223\times {10}^{-6}$

and

Table 6. Table of forcing terms depending on mode number for the tail patch.

	$D_0=D_1\cos \left(\omega t\right)+D_{2}sin\left(\omega t\right)=Dcos\left(\omega t+{\varphi }_1\right)$; where $D=\sqrt{D_1^2+D_2^2}$ ${\mathsf{\Lambda }}_0={\mathsf{\Lambda }}_1{\omega }_p^{2}L\cos \left(\omega t\right)+{\mathsf{\Lambda }}_2{\omega }_p^{2}Lsin\left(\omega t\right)=\mathsf{\Lambda }{\omega }_p^{2}Lcos\left(\omega t+{\varphi }_2\right)$; where $\mathsf{\Lambda }=\sqrt{{\mathsf{\Lambda }}_1^2+{\mathsf{\Lambda }}_2^2}$
	Case A	Case B	Case C	Case D	Case E	Case F
$i=1$	$D=30.3592$ $\mathsf{\Lambda }=9.2471\times {10}^{-6}$	$D=15.4722$ $\mathsf{\Lambda }=2.1167\times {10}^{-5}$	$D=22.1269$ $\mathsf{\Lambda }=1.7064\times {10}^{-5}$	$D=20.8087$ $\mathsf{\Lambda }=1.7440\times {10}^{-4}$	$D=0$ $\mathsf{\Lambda }=1.3018\times {10}^{-4}$	$D=0$ $\mathsf{\Lambda }=6.0188\times {10}^{-5}$
$i=2$	$D=358.4153$ $\mathsf{\Lambda }=8.0113\times {10}^{-5}$	$D=286.2125$ $\mathsf{\Lambda }=7.1840\times {10}^{-5}$	$D=278.1450$ $\mathsf{\Lambda }=2.9552\times {10}^{-5}$	$D=84.7969$ $\mathsf{\Lambda }=1.9829\times {10}^{-5}$	$\mathsf{\Lambda }=3.4405\times {10}^{-6}$	$\mathsf{\Lambda }=3.0193\times {10}^{-5}$
$i=3$	$D=1.5766\times {10}^3$ $\mathsf{\Lambda }=9.5397\times {10}^{-5}$	$D=1.4662\times {10}^3$ $\mathsf{\Lambda }=6.3837\times {10}^{-5}$	$D=850.3386$ $\mathsf{\Lambda }=5.0922\times {10}^{-5}$	$D=170.8961$ $\mathsf{\Lambda }=1.3764\times {10}^{-5}$	$\mathsf{\Lambda }=8.1347\times {10}^{-6}$	$\mathsf{\Lambda }=1.7421\times {10}^{-5}$
$i=4$	$D=4.1452\times {10}^3$ $\mathsf{\Lambda }=2.7765\times {10}^{-5}$	$D=3.1095\times {10}^3$ $\mathsf{\Lambda }=2.9907\times {10}^{-5}$	$D=841.8907$ $\mathsf{\Lambda }=2.2707\times {10}^{-5}$	$D=292.5798$ $\mathsf{\Lambda }=5.3014\times {10}^{-6}$	$\mathsf{\Lambda }=7.7336\times {10}^{-6}$	$\mathsf{\Lambda }=1.2410\times {10}^{-5}$
$i=5$	$D=5.2101\times {10}^3$ $\mathsf{\Lambda }=2.2700\times {10}^{-5}$	$D=2.6451\times {10}^3$ $\mathsf{\Lambda }=8.2872\times {10}^{-6}$	$D=912.0545$ $\mathsf{\Lambda }=1.3641\times {10}^{-5}$	$D=408.3896$ $\mathsf{\Lambda }=7.9469\times {10}^{-6}$	$\mathsf{\Lambda }=6.7128\times {10}^{-6}$	$\mathsf{\Lambda }=9.4904\times {10}^{-6}$
$i=6$	$D=5.7004\times {10}^3$ $\mathsf{\Lambda }=3.6307\times {10}^{-6}$	$D=3.8640\times {10}^3$ $\mathsf{\Lambda }=9.8558\times {10}^{-6}$	$D=1.2206\times {10}^3$ $\mathsf{\Lambda }=1.1583\times {10}^{-5}$	$D=433.8227$ $\mathsf{\Lambda }=5.3117\times {10}^{-6}$	$\mathsf{\Lambda }=5.9012\times {10}^{-6}$	$\mathsf{\Lambda }=7.7706\times {10}^{-6}$
$i=7$	$D=6.9111\times {10}^3$ $\mathsf{\Lambda }=8.8546\times {10}^{-6}$	$D=3.8624\times {10}^3$ $\mathsf{\Lambda }=8.0644\times {10}^{-6}$	$D=1.1560\times {10}^3$ $\mathsf{\Lambda }=8.0925\times {10}^{-6}$	$D=379.4358$ $\mathsf{\Lambda }=5.6015\times {10}^{-6}$	$\mathsf{\Lambda }=5.2992\times {10}^{-6}$	$\mathsf{\Lambda }=6.6733\times {10}^{-6}$

. These terms are dependent on the number of modes, and each of these terms is identified by its mode number (i). These terms are found by separating the prescribed motion components that depend on x and time. Magnitudes of these constant terms from the length evaluation are denoted for the double prime term as D and from the double dot term as $\mathsf{\Lambda }$. For case A prescribed motion, D has the largest magnitude, and the tail configuration offers the largest variation of these values, as shown in Table 6. Cases B–D all have an incrementally reduced magnitude when components of the envelope or K are removed from the prescribed motion for head and tail patches. Typically, $\mathsf{\Lambda }$ is found for base-excited systems and for excitations that are out of resonance; it can be seen for all cases that these magnitudes are small for all values of i. The variation of the $\mathsf{\Lambda }$ magnitude for the head case is smaller than the tail case, as is true for D.

Further investigation was carried out to try and explain the impacts of the terms in Table 5 and Table 6. Equation (16) contains a coupling term to Equation (15) along with another excitation from the prescribed motion. The prescribed motion excitation term is independent of the mode number and is known to be able to give us the uncoupled case if the damping is very high. This coupling term is of interest as it has been shown to impact the harvested power and is influenced by the number of modes, this term is defined as

$$\gamma ={\left({{\displaystyle \int }}_{L_1}^{L_2}\left\{{\displaystyle \sum }_{i=1}^N{\varphi }_i^{″}{\dot{q_{𝚤}}}_{ }\right\}dx\right)}_{ }$$

(22)

The piezoelectric coupling term $\theta$ is not considered for the definition of this term. The effects of case D motion on $\gamma$ are shown in Figure 10. This figure shows the maximum value of these time-dependent terms to determine the $N$ mode number effect for a range of frequencies. Figure 10 demonstrates that there are impacts on the convergence in Equation (16) for the different configurations of the patch under case D motion. Figure 10a shows that the variation in the maximal values of $\gamma$ for the head patch configuration converges with a small variation when increasing the number of modes. Figure 10b demonstrates that, for the tail patch attachment, there is a clear convergence for $N=6$ and 7 with a large variation between the modes under consideration.

Table 7. Values for ${\gamma }_{max}$ for $f_p$ = 5 Hz and $\xi =0.3$.

	Head $\mathit{R}={10}^3 \mathbf{\Omega }$	Head $\mathit{R}={10}^5 \mathbf{\Omega }$	Head $\mathit{R}={10}^8 \mathbf{\Omega }$	Tail $\mathit{R}={10}^3 \mathbf{\Omega }$	Tail $\mathit{R}={10}^5 \mathbf{\Omega }$	Tail $\mathit{R}={10}^8 \mathbf{\Omega }$
Case A
$N=1$	4.977	4.969	$4.775$	0.2392	0.2404	0.1825
$N=2$	5.589	5.574	$5.604$	1.862	1.868	1.523
$N=3$	5.494	5.479	$5.605$	7.071	7.086	7.121
$N=4$	2.648	2.629	$2.826$	14.727	14.740	14.884
$N=5$	4.661	4.639	$4.854$	19.109	19.118	19.349
$N=6$	4.244	4.227	$4.472$	21.213	21.222	21.412
$N=7$	3.769	3.754	$4.030$	21.705	21.725	21.909
Case D
$N=1$	3.401	3.396	$3.271$	0.164	0.163	0.157
$N=2$	2.809	2.811	$2.691$	0.587	0.586	0.555
$N=3$	2.823	2.821	$2.693$	1.175	1.171	1.110
$N=4$	3.037	3.042	$2.889$	1.788	1.785	1.695
$N=5$	2.959	2.962	$2.811$	2.248	2.241	2.130
$N=6$	2.992	2.996	$2.845$	2.432	2.438	2.329
$N=7$	3.031	3.035	$2.881$	2.479	2.477	2.369

gives insight into $\gamma$ for different load resistances and cases of motion for $N=1\to 7$. The unclear convergence phenomenon of the maximal value of $\gamma$ for the head patch can be observed for the case A values for the three load resistances used and in case D for the two other load resistances. The transition of results to convergence in the tail patch case is perceived for the three electrical load resistances for case D. This suggests that seven modes are possibly enough to approximate the performance of the system for case D. Case A could potentially present issues for the power harvester due to the undulatory terms impacting the system and causing a larger variation in the tail patch configuration. The load resistance has been proven to have some impact on the convergence of these systems, but it is considered minimal compared to the case of motion.

5.2. Case Motion Impacts the Interaction between Prescribed and Relative Motions

During the convergence analysis investigation, $v_r\left(L\right)$ for case A motion is perceived to be larger than the maximal prescribed value at this location. During case D motion, $v_r\left(L\right)$ also increased relative to the same value of the prescribed motion. It is necessary to determine how the total displacement is affected to distinguish if the system is moving outside the bounds of the desired envelope of the motion. Time histories are used to help explain this phenomenon, as shown in Figure 11. It is obvious in Figure 11 that the relative displacement is larger for cases A and D for the two patch configurations. However, the phase of $v_r$ results in a more constructive build in case A and destructive loss in case D compared to the prescribed motion. This results in a total amplitude ($v_t$) that is larger than $v_p$ for case A and smaller in case D.

Table 8. Maximum $v\left(L\right) \left(\mathrm{m}\right)$ amplitude for $R={10}^5 \mathsf{\Omega }$ and $\xi =0.3$.

	Head			Tail
	Case A
Frequency (Hz)	Prescribed ($v_p$) $\times {10}^{-2}$	Relative ($v_r$) $\times {10}^{-2}$	Total motion ($v_t$)$\times {10}^{-2}$	Prescribed ($v_p$) $\times {10}^{-2}$	Relative ($v_r$) $\times {10}^{-2}$	Total motion ($v_t$)$\times {10}^{-2}$
1	1	2.1174	1.4586	1	2.1230	1.4580
2	1	2.1232	1.4558	1	2.1365	1.4571
3	1	2.1298	1.4355	1	2.1507	1.4574
4	1	2.1371	1.4510	1	2.1662	1.4603
5	1	2.1453	1.4488	1	2.1824	1.4623
6	1	2.1534	1.4524	1	2.2017	1.4683
7	1	2.1634	1.4528	1	2.2212	1.4749
8	1	2.1704	1.4457	1	2.2347	1.4828
9	1	2.1827	1.4568	1	2.2639	1.4934
	Case D
	Prescribed ($v_p$) $\times {10}^{-2}$	Relative ($v_r$) $\times {10}^{-2}$	Total motion ($v_t$)$\times {10}^{-3}$	Prescribed ($v_p$) $\times {10}^{-2}$	Relative ($v_r$) $\times {10}^{-2}$	Total motion ($v_t$)$\times {10}^{-3}$
1	1	1.62093	6.21132	1	1.62121	6.21545
2	1	1.62096	6.21779	1	1.62121	6.22787
3	1	1.62175	6.22957	1	1.62256	6.24913
4	1	1.62145	6.24256	1	1.62174	6.28299
5	1	1.62169	6.26169	1	1.62220	6.31971
6	1	1.62418	6.28476	1	1.62745	6.36952
7	1	1.62488	6.30981	1	1.62995	6.42618
8	1	1.62523	6.33117	1	1.63273	6.47224
9	1	1.62769	6.38672	1	1.63534	6.56591

details the values for prescribed, relative and total displacement for varying prescribed frequencies. The larger total motion for case A remains true for the frequency range, as well the smaller total motion in case D.

Similar to the displacement time history investigation, an investigation into $\gamma$ and the prescribed term in Equation (16) is carried out to determine how the coupled motion case affects the interaction of the terms in this equation. Figure 12 details the response of the coupled systems for $f_p=5 \mathrm{Hz}$. These time histories are observed for the amplitude and phase for the $\gamma$ and prescribed term. Figure 12a,c demonstrates that the amplitude for $\gamma$ and the prescribed term are very similar for the tail patch configurations. For the tail configurations, the amplitude for these terms vary, as shown in Figure 12b,d. The prescribed term and $\gamma$ tend to be out of phase for the head patch and in phase for the tail patch. The values of the amplitude are detailed in

Table 9. Current investigation in Gauss law terms for $R={10}^5 \mathsf{\Omega }$ and $\xi =0.3$.

	Head			Tail
	Case A
Frequency (Hz)	${\left\{{\left.{\dot{v_p}}^{\prime }\right|}_0^{L_1}\right\}}_{RMS}$	${\left\{\gamma \right\}}_{RMS}$	Total $i_c{}_{RMS}$	${\left\{{\left.{\dot{v_p}}^{\prime }\right|}_{L_1}^L\right\}}_{RMS}$	${\left\{\gamma \right\}}_{RMS}$	Total $i_c{}_{RMS}$
1	0.8288	0.7578	0.0726	5.3298	2.9888	7.7665
2	1.6575	1.5163	0.1508	10.6585	6.0014	15.4966
3	2.4865	2.2633	0.2252	15.9895	9.0109	23.2152
4	3.3151	3.0114	0.3174	21.3006	12.0326	30.9158
5	4.1441	3.7537	0.4406	26.6258	15.0815	38.5161
6	4.9727	4.4988	0.6031	31.9756	18.1303	46.2722
7	5.8018	5.2470	0.7164	37.3089	21.2052	53.9468
8	6.6302	5.9799	0.8135	42.6012	24.3132	61.5902
9	7.4595	6.7297	1.0927	47.9686	27.4530	69.1966
	Case D
	${\left\{{\left.{\dot{v_p}}^{\prime }\right|}_0^{L_1}\right\}}_{RMS}$	${\left\{\gamma \right\}}_{RMS}$	Total $i_c{}_{RMS}$	${\left\{{\left.{\dot{v_p}}^{\prime }\right|}_{L_1}^L\right\}}_{RMS}$	${\left\{\gamma \right\}}_{RMS}$	Total $i_c{}_{RMS}$
1	0.6126	0.6049	0.0153	0.6126	1.4386	0.8263
2	1.2252	1.2106	0.0417	1.2252	2.8771	1.6524
3	1.8378	1.8199	0.0764	1.8378	4.3324	2.4954
4	2.4456	2.4317	0.1306	2.4456	5.7776	3.3320
5	3.0631	3.0347	0.1959	3.0630	7.1980	4.1493
6	3.6757	3.6473	0.2639	3.6756	8.6734	5.0050
7	4.2882	4.2407	0.3434	4.2882	10.1217	5.8525
8	4.8912	4.8481	0.4551	4.8912	11.5482	6.7341
9	5.5134	5.4692	0.5765	5.5134	13.0296	7.5890

for varying frequencies.

5.3. Impacts of Load Resistance and Damping on the System’s Performance

All previous studies were done using a fixed load resistance or damping ratio. The two configurations of patches and motion are used to investigate their performance impacts on the BCFEH when responding to a variation of these terms, as shown in Figure 13. A performance table is then presented to detail characteristics of the respective configuration. The performance of the BCFEH is known to be highly dependent on the kind of motion. Figure 13a,c,e shows the $v_r\left(L\right)$ response to the variation of the load resistance. The relative motion undergoes the shunt damping phenomenon when the load resistance for optimal power harvesting occurs. Figure 13b,d,f presents the harvested power for the four configurations. The relative motion for the head patch in case A shows an increase in amplitude after the shunt damping impact, which can be attributed to the undulatory term affecting the motion. The tail patch in case A demonstrates a reduction in amplitude and does not exhibit the same phenomenon; thus, the structural stiffness over the length possibly has an impact. The relative motion in case D is more impacted by different $\xi$ values than the power harvested, as shown in Figure 13e,f, respectively. The head and tail patches undergoing case D motion are very similar, and so only the tail patch results are shown. Values for the case A motion can be found in

Table 10. Case A motion values of $v_r\left(L\right)$ and $P_{avg}$ with corresponding resistances for different $\xi$ and patch configuration when $f=5 \mathrm{Hz}$.

Damping Value	Case A Head		Case A Tail
$\xi =0.1$	$v_r{\left(L\right)}_{ }{}_{max}=2.1354\times {10}^{-2} m$ $v_r{\left(L\right)}_{ }{}_{min}=1.9605\times {10}^{-2} m$ ${P_{avg}}_{max}=3.5709\times {10}^{-3} W$ ${P_{avg}}_{Unc}=1.6347\times {10}^{-3} W$	$R_{{\mathrm{v}}_{\mathrm{r}}{}_{max}}=5.1506\times {10}^2 \mathsf{\Omega }$ $R_{{\mathrm{v}}_{\mathrm{r}}{}_{min}}=2.3598\times {10}^6 \mathsf{\Omega }$ $R_{P_{max}}=2.3598\times {10}^6 \mathsf{\Omega }$ ${R_P}_{Unc}=2.3598\times {10}^6 \mathsf{\Omega }$	$v_r{\left(L\right)}_{ }{}_{max}=2.1539\times {10}^{-2} m$ $v_r{\left(L\right)}_{ }{}_{min}=1.7670\times {10}^{-2} m$ ${P_{avg}}_{max}=4.3121\times {10}^{-2} W$ ${P_{avg}}_{Unc}=6.7598\times {10}^{-2} W$	$R_{{\mathrm{v}}_{\mathrm{r}}{}_{max}}=2.8683\times {10}^5 \mathsf{\Omega }$ $R_{{\mathrm{v}}_{\mathrm{r}}{}_{min}}=4.9535\times {10}^7 \mathsf{\Omega }$ ${R_P}_{max}=2.3598\times {10}^6 \mathsf{\Omega }$ ${R_P}_{Unc}=2.3598\times {10}^6 \mathsf{\Omega }$
$\xi =0.3$	$v_r{\left(L\right)}_{ }{}_{max}=2.1596\times {10}^{-2} m$ $v_r{\left(L\right)}_{ }{}_{min}=1.9809\times {10}^{-2} m$ ${P_{avg}}_{max}=3.6815\times {10}^{-3} W$ ${P_{avg}}_{Unc}=1.6347\times {10}^{-3} W$	$R_{{\mathrm{v}}_{\mathrm{r}}{}_{max}}=4.0753\times {10}^2 \mathsf{\Omega }$ $R_{{\mathrm{v}}_{\mathrm{r}}{}_{min}}=2.3598\times {10}^6 \mathsf{\Omega }$ ${R_P}_{max}=2.3598\times {10}^6 \mathsf{\Omega }$ $R_P{}_{Unc}=2.3598\times {10}^6 \mathsf{\Omega }$	$v_r{\left(L\right)}_{ }{}_{max}=2.1821\times {10}^{-2} m$ $v_r{\left(L\right)}_{ }{}_{min}=1.7929\times {10}^{-2} m$ ${P_{avg}}_{max}=4.3041\times {10}^{-2} W$ ${P_{avg}}_{Unc}=6.7598\times {10}^{-2} W$	$R_{{\mathrm{v}}_{\mathrm{r}}{}_{max}}=1.4208\times {10}^5 \mathsf{\Omega }$ $R_{{\mathrm{v}}_{\mathrm{r}}{}_{min}}=1\times {10}^9 \mathsf{\Omega }$ ${R_P}_{max}=2.3598\times {10}^6 \mathsf{\Omega }$ ${R_P}_{Unc}=2.3598\times {10}^6 \mathsf{\Omega }$
$\xi =0.5$	$v_r{\left(L\right)}_{ }{}_{max}=2.1795\times {10}^{-2} m$ $v_r{\left(L\right)}_{ }{}_{min}=2.0035\times {10}^{-2} m$ ${P_{avg}}_{max}=3.7978\times {10}^{-3} W$ ${P_{avg}}_{Unc}=1.6347\times {10}^{-3} W$	$R_{{\mathrm{v}}_{\mathrm{r}}{}_{max}}=5.1506\times {10}^2 \mathsf{\Omega }$ $R_{{\mathrm{v}}_{\mathrm{r}}{}_{min}}=2.3598\times {10}^6 \mathsf{\Omega }$ $R_{P_{max}}=2.3598\times {10}^6 \mathsf{\Omega }$ ${R_P}_{Unc}=2.3598\times {10}^6 \mathsf{\Omega }$	$v_r{\left(L\right)}_{ }{}_{max}=2.2057\times {10}^{-2} m$ $v_r{\left(L\right)}_{ }{}_{min}=1.8238\times {10}^{-2} m$ ${P_{avg}}_{max}=4.2499\times {10}^{-2} W$ ${P_{avg}}_{Unc}=6.7598\times {10}^{-2} W$	$R_{{\mathrm{v}}_{\mathrm{r}}{}_{max}}=2.8683\times {10}^5 \mathsf{\Omega }$ $R_{{\mathrm{v}}_{\mathrm{r}}{}_{min}}=1\times {10}^9 \mathsf{\Omega }$ ${R_P}_{max}=2.3598\times {10}^6 \mathsf{\Omega }$ $R_P{}_{Unc}=2.3598\times {10}^6 \mathsf{\Omega }$

, and values for case D motion are shown in

Table 11. Case D values of $v_r\left(L\right)$ and $P_{avg}$ with corresponding resistances for different $\xi$ and patch configuration when $f=5 \mathrm{Hz}$.

Damping Value	Case D Head		Case D Tail
$\xi =0.1$	$v_r{\left(L\right)}_{ }{}_{min}=1.6262\times {10}^{-2} m$ $v_r{\left(L\right)}_{ }{}_{min}=1.5277\times {10}^{-2} m$ ${P_{avg}}_{max}=3.7462\times {10}^{-3} W$ ${P_{avg}}_{Unc}=8.9291\times {10}^{-4} W$	$R_{{\mathrm{v}}_{\mathrm{r}}{}_{max}}=1.6103\times {10}^4 \mathsf{\Omega }$ $R_{{\mathrm{v}}_{\mathrm{r}}{}_{min}}=4.5204\times {10}^7\mathsf{\Omega }$ $R_{P_{max}}=2.3598\times {10}^6 \mathsf{\Omega }$ $R_P{}_{Unc}= 2.3598\times {10}^6 \mathsf{\Omega }$	$v_r{\left(L\right)}_{ }{}_{max}=2.1823\times {10}^{-2} m$ $v_r{\left(L\right)}_{ }{}_{min}=1.7929\times {10}^{-2} m$ ${P_{avg}}_{max}=3.7788\times {10}^{-3} W$ ${P_{avg}}_{Unc}=8.9291\times {10}^{-4} W$	$R_{{\mathrm{v}}_{\mathrm{r}}{}_{max}}=1.2690\times {10}^5 \mathsf{\Omega }$ $R_{{\mathrm{v}}_{\mathrm{r}}{}_{min}}=1\times {10}^8 \mathsf{\Omega }$ $R_{P_{max}}=2.3598\times {10}^6 \mathsf{\Omega }$ $R_P{}_{Unc}=2.3598\times {10}^6 \mathsf{\Omega }$
$\xi =0.2$	$v_r{\left(L\right)}_{ }{}_{max}=1.6246\times {10}^{-2} m$ $v_r{{\left(L\right)}_{avg}}_{min}=1.5266\times {10}^{-2} m$ ${P_{avg}}_{max}=3.7588\times {10}^{-3}W$ ${P_{avg}}_{Unc}=8.9291\times {10}^{-4} W$	$R_{{\mathrm{v}}_{\mathrm{r}}{}_{max}}=6.7233\times {10}^2 \mathsf{\Omega }$ $R_{{\mathrm{v}}_{\mathrm{r}}{}_{min}}=6.2102\times {10}^7\mathsf{\Omega }$ $R_P{}_{max}=2.3598\times {10}^6 \mathsf{\Omega }$ $R_P{}_{Unc}=2.3598\times {10}^6 \mathsf{\Omega }$	$v_r{\left(L\right)}_{ }{}_{max}=1.6276\times {10}^{-2} m$ $v_r{\left(L\right)}_{ }{}_{min}=1.5106\times {10}^{-2} m$ ${P_{avg}}_{max}=3.7617\times {10}^{-3} W$ ${P_{avg}}_{Unc}=8.9291\times {10}^{-4} W$	$R_{{\mathrm{v}}_{\mathrm{r}}{}_{max}}=1.0826\times {10}^3$ $R_{{\mathrm{v}}_{\mathrm{r}}{}_{min}}=3.8566\times {10}^7 \mathsf{\Omega }$ $R_{P_{max}}=2.3598\times {10}^6 \mathsf{\Omega }$ $R_P{}_{Unc}=2.3598\times {10}^6 \mathsf{\Omega }$
$\xi =0.3$	$v_r{\left(L\right)}_{ }{}_{max}=1.6222\times {10}^{-2} m$ $v_r{\left(L\right)}_{ }{}_{min}=1.5248\times {10}^{-2} m$ ${P_{avg}}_{max}=3.7245\times {10}^{-3} W$ ${P_{avg}}_{Unc}=8.9291\times {10}^{-4} W$	$R_{{\mathrm{v}}_{\mathrm{r}}{}_{max}}=100 \mathsf{\Omega }$ $R_{{\mathrm{v}}_{\mathrm{r}}{}_{min}}=3.8566\times {10}^7 \mathsf{\Omega }$ $R_{P_{max}}=2.3598\times {10}^6 \mathsf{\Omega }$ $R_P{}_{Unc}=2.3598\times {10}^6 \mathsf{\Omega }$	$v_r{\left(L\right)}_{ }{}_{max}=1.6233\times {10}^{-2} m$ $v_r{\left(L\right)}_{ }{}_{min}=1.5070\times {10}^{-2} m$ ${P_{avg}}_{max}=3.7181\times {10}^{-3} W$ ${P_{avg}}_{Unc}=8.9291\times {10}^{-4} W$	$R_{{\mathrm{v}}_{\mathrm{r}}{}_{max}}=1.6102\times {10}^2 \mathsf{\Omega }$ $R_{{\mathrm{v}}_{\mathrm{r}}{}_{min}}=3.8566\times {10}^7\mathsf{\Omega }$ $R_{P_{max}}=2.3598\times {10}^6 \mathsf{\Omega }$ $R_P{}_{Unc}=2.3598\times {10}^6 \mathsf{\Omega }$
$\xi =0.4$	$v_r{\left(L\right)}_{ }{}_{max}=1.6194\times {10}^{-2} m$ $v_r{\left(L\right)}_{ }{}_{min}=1.5216\times {10}^{-2} m$ ${P_{avg}}_{max}=3.7091\times {10}^{-3} W$ ${P_{avg}}_{Unc}=8.9291\times {10}^{-3} W$	$R_{{\mathrm{v}}_{\mathrm{r}}{}_{max}}=2.8072\times {10}^3 \mathsf{\Omega }$ $R_{{\mathrm{v}}_{\mathrm{r}}{}_{min}}=6.2102\times {10}^7 \mathsf{\Omega }$ $R_{P_{max}}=2.3598\times {10}^6 \mathsf{\Omega }$ $R_P{}_{Unc}=2.3598\times {10}^6 \mathsf{\Omega }$	$v_r{\left(L\right)}_{ }{}_{max}=1.6197\times {10}^{-2} m$ $v_r{\left(L\right)}_{ }{}_{min}=1.5034\times {10}^{-2} m$ ${P_{avg}}_{max}=3.7117\times {10}^{-3} W$ ${P_{avg}}_{Unc}=8.9291\times {10}^{-4} W$	$R_{{\mathrm{v}}_{\mathrm{r}}{}_{max}}=1.0826\times {10}^3 \mathsf{\Omega }$ $R_{{\mathrm{v}}_{\mathrm{r}}{}_{min}}=6.2102\times {10}^7 \mathsf{\Omega }$ $R_{P_{max}}=2.3598\times {10}^6 \mathsf{\Omega }$ $R_P{}_{Unc}=2.3598\times {10}^6 \mathsf{\Omega }$
$\xi =0.5$	$v_r{\left(L\right)}_{ }{}_{max}=1.6156\times {10}^{-2} m$ $v_r{\left(L\right)}_{ }{}_{min}=1.5191\times {10}^{-2} m$ ${P_{avg}}_{max}=3.7048\times {10}^{-3} W$ ${P_{avg}}_{Unc}=8.9291\times {10}^{-4} W$	$R_{{\mathrm{v}}_{\mathrm{r}}{}_{max}}=1.7433\times {10}^3 \mathsf{\Omega }$ $R_{{\mathrm{v}}_{\mathrm{r}}{}_{min}}=2.3950\times {10}^7 \mathsf{\Omega }$ $R_{P_{max}}=2.3598\times {10}^6 \mathsf{\Omega }$ $R_P{}_{Unc}=2.3598\times {10}^6 \mathsf{\Omega }$	$v_r{\left(L\right)}_{ }{}_{max}=1.6176\times {10}^{-2} m$ $v_r{\left(L\right)}_{ }{}_{min}=1.5020\times {10}^{-2} m$ ${P_{avg}}_{max}=3.6945\times {10}^{-3} W$ ${P_{avg}}_{Unc}=8.9291\times {10}^{-4} W$	$R_{{\mathrm{v}}_{\mathrm{r}}{}_{max}}=2.8072\times {10}^3 \mathsf{\Omega }$ $R_{{\mathrm{v}}_{\mathrm{r}}{}_{min}}=1\times {10}^8 \mathsf{\Omega }$ $R_{P_{max}}=2.3598\times {10}^6 \mathsf{\Omega }$ $R_P{}_{Unc}=2.3598\times {10}^6 \mathsf{\Omega }$
$\xi =0.6$	$v_r{\left(L\right)}_{ }{}_{max}=1.6148\times {10}^{-2} m$ $v_r{\left(L\right)}_{ }{}_{min}=1.5179\times {10}^{-2} m$ ${P_{avg}}_{max}=3.6970\times {10}^{-3} W$ ${P_{avg}}_{Unc}=8.9291\times {10}^{-4} W$	$R_{{\mathrm{v}}_{\mathrm{r}}{}_{max}}=2.5929\times {10}^2 \mathsf{\Omega }$ $R_{{\mathrm{v}}_{\mathrm{r}}{}_{min}}=1\times {10}^8 \mathsf{\Omega }$ $R_{P_{max}}=2.3598\times {10}^6 \mathsf{\Omega }$ $R_P{}_{Unc}=2.3598\times {10}^6 \mathsf{\Omega }$	$v_r{\left(L\right)}_{ }{}_{max}=1.6138\times {10}^{-2} m$ $v_r{\left(L\right)}_{ }{}_{min}=1.4985\times {10}^{-2} m$ ${P_{avg}}_{max}=3.6741\times {10}^{-3} W$ ${P_{avg}}_{Unc}=8.9291\times {10}^{-4} W$	$R_{{\mathrm{v}}_{\mathrm{r}}{}_{max}}=4.1753\times {10}^2 \mathsf{\Omega }$ $R_{{\mathrm{v}}_{\mathrm{r}}{}_{min}}=3.8566\times {10}^7 \mathsf{\Omega }$ $R_{P_{max}}=2.3598\times {10}^6 \mathsf{\Omega }$ $R_P{}_{Unc}=2.3598\times {10}^6 \mathsf{\Omega }$
$\xi =0.7$	$v_r{\left(L\right)}_{ }{}_{max}=1.6113\times {10}^{-2} m$ $v_r{\left(L\right)}_{ }{}_{min}=1.5166\times {10}^{-2} m$ ${P_{avg}}_{max}=3.6842\times {10}^{-3} W$ ${P_{avg}}_{Unc}=8.9291\times {10}^{-4} W$	$R_{{\mathrm{v}}_{\mathrm{r}}{}_{max}}=2.8072\times {10}^3 \mathsf{\Omega }$ $R_{{\mathrm{v}}_{\mathrm{r}}{}_{min }}=2.3950\times {10}^7 \mathsf{\Omega }$ $R_{P_{max}}=2.3598\times {10}^6 \mathsf{\Omega }$ $R_P{}_{Unc}=2.3598\times {10}^6 \mathsf{\Omega }$	$v_r{\left(L\right)}_{ }{}_{max}=1.6084\times {10}^{-2} m$ $v_r{\left(L\right)}_{ }{}_{min}=1.4944\times {10}^{-2} m$ ${P_{avg}}_{max}=3.6434\times {10}^{-3} W$ ${P_{avg}}_{Unc}=8.9291\times {10}^{-4} W$	$R_{{\mathrm{v}}_{\mathrm{r}}{}_{max}}=1.6102\times {10}^2 \mathsf{\Omega }$ $R_{{\mathrm{v}}_{\mathrm{r}}{}_{min}}=2.3950\times {10}^6 \mathsf{\Omega }$ $R_{P_{max}}=2.3598\times {10}^6 \mathsf{\Omega }$ $R_P{}_{Unc}=2.3598\times {10}^6 \mathsf{\Omega }$

.

Case A motion presents harvested power curves which should be considered to be uncertain, especially for the tail patch configuration. The power harvested for the head patch in case A is dramatically higher than the uncoupled power harvester results, which can be seen in the convergence investigations. Compared to the tail patch in case A, the head patch could potentially harvest a similar amount to the tail patch configuration, as shown in Table 10. The harvested power appears to be very similar for both patch configurations case D motion, as indicated in Table 11.

The importance of damping on the relative displacement shown in Table 10 and Table 11 needs some more detailed discussion. The head patch placement for the BCFEH under case D motion considering $\xi =0.1$ shows a maximum relative displacement amplitude of $1.6262\times {10}^{-2} m$. When $\xi =0.7$, this maximum amplitude is $1.6113\times {10}^{-2} m$. It is clear that the increase in damping results in a decrease in amplitude. Therefore, the minimum relative displacement values for $\xi =0.1$ and 0.7 are $1.5277\times {10}^{-2} m$ and $1.5167\times {10}^{-2} m$, respectively, and thus follow the same logic. Regarding the tail patch placement for the BFCEH under case D motion, the maximum relative displacement is more affected by the change of damping. For values of $\xi$ equal to 0.1 and 0.7, relative displacement values of $2.1823\times {10}^{-2} m$ and $1.6084\times {10}^{-2} m$, respectively, are obtained. The minimum relative displacements of the tail tip are $1.7929\times {10}^{-2} m$ and $1.4944\times {10}^{-2} m$. The tail patch configuration is more sensitive to the variations of the damping that affects the tip displacement.

The harvested power is only slightly affected by the damping variation. The kind of motion and patch configuration do not largely affect the load resistance response of the energy harvester. The slight variation can be seen for the average powers for the head patch and tail patch, considering $\xi =0.1$ and 0.7, which are $3.7462\times {10}^{-3} W$ and $3.6842\times {10}^{-3} W$ and $3.7788\times {10}^{-3} W$ and $3.6434\times {10}^{-3} W$, respectively. Since the uncoupled system is unaffected by damping and the energy harvester exhibits a very similar interaction with motion, the maximum $P_{avg}$ for the uncoupled head and tail patch configuration is $8.9291\times {10}^{-4} W$. Comparing the coupled and uncoupled powers, the case D motion patch configurations are shown to be able to harvest more power. The head patch coupled case exhibits less variation from the uncoupled power. The tail patch configuration is more sensitive to the damping and presents cases where it harvests less power for higher values of damping.

Investigating the load resistances which correspond to these values, it is clear that there is an equal optimal load resistance for the coupled and uncoupled harvested power. The resistance related to the maximum or minimum relative displacement does not have this unity. Previous energy harvester works have shown that the resistance that relates to maximal power is an optimal resistance of $R_{opt}\approx 1/C_p\omega$. For the two patch configurations forced with $f_p=5 \mathrm{Hz}$, $R_{opt}\approx 2.3959\times {10}^6 \mathsf{\Omega }$ for both cases as the same length of patch is used in both configurations, and thus $C_p$ is equivalent.

To investigate how the change in $\xi$ could possibly alter the total motion, the methods employed in Section 5.3 are reutilized for three values of damping. The results are organized into

Table 12. Impacts on max $v\left(L\right)$ for varying $\xi$ when $R={10}^5 \mathsf{\Omega }$ and $f_p=5 \mathrm{Hz}$.

	Head			Tail
	Case A
$\xi$	Prescribed ($v_p$) $\times {10}^{-2}$	Relative ($v_r$) $\times {10}^{-2}$	Total motion ($v_t$)$\times {10}^{-2}$	Prescribed ($v_p$) $\times {10}^{-2}$	Relative ($v_r$) $\times {10}^{-2}$	Total motion ($v_t$)$\times {10}^{-2}$
0.1	1	2.1207	1.4587	1	2.1525	1.4736
0.3	1	2.1453	1.4488	1	2.1824	1.4623
0.5	1	2.1654	1.4437	1	2.2203	1.4488
	Case D
	Prescribed ($v_p$) $\times {10}^{-2}$	Relative ($v_r$) $\times {10}^{-2}$	Total motion ($v_t$)$\times {10}^{-3}$	Prescribed ($v_p$) $\times {10}^{-2}$	Relative ($v_r$) $\times {10}^{-2}$	Total motion ($v_t$)$\times {10}^{-3}$
0.1	1	1.6251	6.2513	1	1.6295	6.2953
0.3	1	1.6217	6.2617	1	1.6222	6.3197
0.5	1	1.6150	6.2735	1	1.6165	6.3477

. The impact of the damping on the total motion is not dramatic.

5.4. Prescribed Frequency Impacts on Harvester Performance

It is known that frequency can have an impact on the performance of the system when the BCFEH is operating at prescribed excitations. In this section, we aim to explain how different operating conditions could alter the displacement and power along with identifying load resistances which relate to these values. Figure 14 organizes this investigation for the combinations of patch placement and motion, as done for previous investigations. Additionally, uncoupled harvested power results are compared to the coupled results as in Section 5.3. Again, the head and tail patches undergoing case D motion exhibit very similar results, and so only the tail patch case is shown.

Indeed, the phenomena seen for $f_p=5 \mathrm{Hz}$ from the damping investigation are still present. Two outlier prescribed frequency values are used: $f_p=2 \mathrm{Hz}$ for a low operating frequency and $f_p=8 \mathrm{Hz}$ for excessive frequency operation. These three values are used to give a range of responses. Figure 14 shows that the frequency is influential in the values for the relative displacement and power. Relative displacement is shown in Figure 14a,c for the head and tail configuration in case A, and similarly for case D in Figure 14e,f. Power is categorized in the same manner, with values for case A presented in Figure 14b,d and for case D presented in Figure 14f for the head and tail patch, respectively. There is variation in the relative displacement and power for these three prescribed frequencies, which can be more easily perceived in

Table 13. Case A: table of values for $v_r\left(L\right)$ and $P_{avg}$ considering different $f_p$.

	Case A Head		Case A Tail
	$v_r\left(L\right) \left(m\right)$	Power (W)	$v_r\left(L\right)\times {10}^{-2}\left(m\right)$	Power (W)
$f_p=2 \mathrm{Hz}$	(Max) Amp: $2.1288\times {10}^{-2}$ Res: $4.0753\times {10}^2 \mathsf{\Omega }$ (Min) Amp: $1.9577\times {10}^{-2}$ Res: $6.0209\times {10}^6 \mathsf{\Omega }$	(Max) Amp: $1.4356\times {10}^{-2}$ Res: $6.0209\times {10}^6 \mathsf{\Omega }$ (Uncoupled Max) Amp: $6.5269\times {10}^{-4}$ Res: $6.0209\times {10}^6 \mathsf{\Omega }$	(Max) Amp: $2.1390\times {10}^{-2}$ Res: $5.7904\times {10}^5 \mathsf{\Omega }$ (Min) Amp: $1.7606\times {10}^{-2}$ Res: $1\times {10}^8 \mathsf{\Omega }$	(Max) Amp: $1.7226\times {10}^{-2}$ Res: $6.0209\times {10}^6 \mathsf{\Omega }$ (Uncoupled Max) Amp: $2.6988\times {10}^{-2}$ Res: $6.0209\times {10}^6 \mathsf{\Omega }$
$f_p=5 \mathrm{Hz}$	(Max) Amp: $2.1594\times {10}^{-2}$ Res: $4.0753\times {10}^2 \mathsf{\Omega }$ (Min) Amp: $1.9855\times {10}^{-2}$ Res: $2.3598\times {10}^6 \mathsf{\Omega }$	(Max) Amp: $3.6815\times {10}^{-2}$ Res: $2.3598\times {10}^6 \mathsf{\Omega }$ (Uncoupled Max) Amp: $1.5739\times {10}^{-3}$ Res: $2.3595\times {10}^6 \mathsf{\Omega }$	(Max) Amp: $2.1824\times {10}^{-2}$ Res: $1.4208\times {10}^5 \mathsf{\Omega }$ (Min) Amp: $1.7929\times {10}^{-2}$ Res: $1\times {10}^8 \mathsf{\Omega }$	(Max) Amp: $4.2388\times {10}^{-2}$ Res: $2.3598\times {10}^6 \mathsf{\Omega }$ (Uncoupled Max) Amp: $6.5084\times {10}^{-2}$ Res: $2.3598\times {10}^6 \mathsf{\Omega }$
$f_p=8 \mathrm{Hz}$	(Max) Amp: $2.1945\times {10}^{-2}$ Res: $2.5514\times {10}^2 \mathsf{\Omega }$ (Min) Amp: $2.0168\times {10}^{-2}$ Res: $1.4774\times {10}^6 \mathsf{\Omega }$	(Max) Amp: $6.0035\times {10}^{-2}$ Res: $1.4474\times {10}^6 \mathsf{\Omega }$ (Uncoupled Max) Amp: $2.5777\times {10}^{-3}$ Res: $1.4774\times {10}^6 \mathsf{\Omega }$	(Max) Amp: $2.2356\times {10}^{-2}$ Res: $2.7585\times {10}^4 \mathsf{\Omega }$ (Min) Amp: $1.8371\times {10}^{-2}$ Res: $5.6234\times {10}^7 \mathsf{\Omega }$	(Max) Amp: $6.6487\times {10}^{-2}$ Res: $1.4474\times {10}^6 \mathsf{\Omega }$ (Uncoupled Max) Amp: $1.0659\times {10}^{-1}$ Res: $1.4774\times {10}^6 \mathsf{\Omega }$

and

Table 14. Case D: table of values for $v_r\left(L\right)$ and $P_{avg}$ considering different $f_p$.

	Case D Head		Case D Tail
	$v_r\left(L\right) \left(m\right)$	Power (W)	$v_r\left(L\right)\times {10}^{-2}\left(m\right)$	Power (W)
$f_p=2 \mathrm{Hz}$	(Max) Amp: $1.6212\times {10}^{-2}$ Res: $1\times {10}^2 \mathsf{\Omega }$ (Min) Amp: $1.5245\times {10}^{-2}$ Res: $7.9123\times {10}^7 \mathsf{\Omega }$	(Max) Amp:$1.4979\times {10}^{-3}$ Res: $6.0209\times {10}^6 \mathsf{\Omega }$ (Uncoupled Max) Amp: $3.5650\times {10}^{-4}$ Res: $6.0209\times {10}^6 \mathsf{\Omega }$	(Max) Amp: $1.6216\times {10}^{-2}$ Res: $1.1242\times {10}^5 \mathsf{\Omega }$ (Min) Amp: $1.5065\times {10}^{-2}$ Res: $1\times {10}^8 \mathsf{\Omega }$	(Max) Amp:$1.4897\times {10}^{-3}$ Res: $6.0209\times {10}^6 \mathsf{\Omega }$ (Uncoupled Max) Amp: $3.5650\times {10}^{-4}$ Res: $6.0209\times {10}^6 \mathsf{\Omega }$
$f_p=5 \mathrm{Hz}$	(Max) Amp: $1.6225\times {10}^{-2}$ Res: $4.0753\times {10}^2 \mathsf{\Omega }$ (Min) Amp:$1.5247\times {10}^{-2}$ Res: $3.9194\times {10}^7 \mathsf{\Omega }$	(Max) Amp: $3.6557\times {10}^{-3}$ Res: $2.3598\times {10}^6 \mathsf{\Omega }$ (Uncoupled Max) Amp: $8.5969\times {10}^{-4}$ Res: $2.3598\times {10}^6 \mathsf{\Omega }$	(Max) Amp: $1.6325\times {10}^{-2}$ Res: $2.0187\times {10}^2 \mathsf{\Omega }$ (Min) Amp: $1.5069\times {10}^{-2}$ Res: $1.9415\times {10}^7 \mathsf{\Omega }$	(Max) Amp: $3.6695\times {10}^{-3}$ Res: $2.3598\times {10}^6 \mathsf{\Omega }$ (Uncoupled Max) Amp: $8.5969\times {10}^{-4}$ Res: $2.3598\times {10}^6 \mathsf{\Omega }$
$f_p=8 \mathrm{Hz}$	(Max) Amp: $1.6260\times {10}^{-2}$ Res: $3.3529\times {10}^3 \mathsf{\Omega }$ (Min) Amp: $1.5271\times {10}^{-2}$ Res: $3.9194\times {10}^7 \mathsf{\Omega }$	(Max) Amp: $5.7958\times {10}^{-3}$ Res: $1.4774\times {10}^6 \mathsf{\Omega }$ (Uncoupled Max) Amp: $1.4079\times {10}^{-3}$ Res: $1.4774\times {10}^6 \mathsf{\Omega }$	(Max) Amp: $1.6347\times {10}^{-2}$ Res: $5.1506\times {10}^2 \mathsf{\Omega }$ (Min) Amp: $1.5149\times {10}^{-2}$ Res: $4.935\times {10}^7 \mathsf{\Omega }$	(Max) Amp: $5.7732\times {10}^{-3}$ Res: $1.4774\times {10}^6 \mathsf{\Omega }$ (Uncoupled Max) Amp: $1.4079\times {10}^{-3}$ Res: $1.4774\times {10}^6 \mathsf{\Omega }$

for case A and case D motion, respectively.

The impact of the prescribed frequency is readily apparent when looking at the percentage of value variation when $f_p=2$ and 8 Hz. The relative displacement variation is taken between the determined maximal values, while the maximum power is compared for the coupled and uncoupled methods for each of the different frequencies. Case A’s maximum relative displacement percentage variations are 2.99% and 4.32% for the head and tail patch; the minimum value variation is 2.93% and 4.16%. Case D’s variations for the maximum values are 0.3% and 0.8%, and the minimum variations are 0.17% and 0.55% for the head and tail patch configuration at the two outlier frequencies.

The effect on the harvested power is investigated by considering the variation of uncoupled versus coupled values for individual prescribed frequencies, and then uncoupled maximum comparisons for $f_p=2$ and 8 Hz. Case A motion causes the harvested power to be dramatically different to the uncoupled case for the head patch and the power for the tail to be uncertain due to lack of convergence; thus, the percentage variation is not investigated. For case D motion and $f_p=2 \mathrm{Hz}$, the variations of harvested power for the head and tail configuration are 320.17% and 317.87%, values for $f_p=5 \mathrm{Hz}$ are 325.23% and 326.84%, and values for $f_p=8 \mathrm{Hz}$ are 309.05% and 310.06%. The variations of the uncoupled values for $f_p$ values of 2–5 and 5–8 Hz are 59.3% and 36.93% for the head coupled system and 59.4% and 36.44% for the tail coupled system.

As seen in Section 5.3, the load resistance which corresponds to the value for the maximal power in the coupled and uncoupled cases is in unity depending on the prescribed frequency. This continually supports the assumption that $R_{opt}\approx 1/C_p\omega$, which can be considered for the BCFEH regardless of the motion case but is dependent on the patch characteristics and the prescribed frequency. Clearly, an increase in the prescribed frequency is followed by a decrease in the optimal value of the load resistance.

5.5. Prescribed Frequency and Damping with Optimal Resistance

Section 5.3 and Section 5.4 have shown that the optimal resistance can be approximated by $1/C_p{\omega }_p$. Indeed, this value influences the system, but by considering this optimal case, the sensitivity of the BCFEH to the damping and prescribed frequency can be more thoroughly investigated for harvesting under optimal circuitry conditions. Figure 15 shows how these terms vary the performance of the head and tail patch configurations for Case D. Case A motion results for this optimal resistance configuration are not included in Figure 15 but are presented in

Table 15. Maximum and minimum values for case A and D motion and head and tail patch placement for optimal resistances and variable $\xi$ and $f_p$.

	Case A Head	Case A Tail
$v_r\left(L\right) \left(m\right)$	(Max) $2.0601\times {10}^{-2}$ (Min) $1.9444\times {10}^{-2}$	(Max) $2.1044\times {10}^{-2}$ (Min) $1.9549\times {10}^{-2}$
$P_{avg} \left(W\right)$	(Max) $6.5542\times {10}^{-2}$ (Min) $3.5219\times {10}^{-3}$ (Uncoupled Max) $2.6159\times {10}^{-3}$ (Uncoupled Min) $1.6350\times {10}^{-4}$	(Max) $6.9594\times {10}^{-2}$ (Min) $4.2724\times {10}^{-3}$ (Uncoupled Max) $1.0817\times {10}^{-1}$ (Uncoupled Min) $6.7606\times {10}^{-3}$
	Case D Head	Case D Tail
$v_r\left(L\right) \left(m\right)$	(Max) $1.5764\times {10}^{-2}$ (Min) $1.5413\times {10}^{-2}$	(Max) $1.5755\times {10}^{-2}$ (Min) $1.5055\times {10}^{-2}$
$P_{avg} \left(W\right)$	(Max) $6.0246\times {10}^{-3}$ (Min) $3.7043\times {10}^{-4}$ (Uncoupled Max) $1.4288\times {10}^{-3}$ (Uncoupled Min) $8.9302\times {10}^{-5}$	(Max) $6.0971\times {10}^{-3}$ (Min) $3.7046\times {10}^{-4}$ (Uncoupled Max) $1.4288\times {10}^{-3}$ (Uncoupled Min) $8.9302\times {10}^{-5}$

. The relative motion is affected by $\xi$ and $f_p$, but as shown in Section 5.3 and Section 5.4, these values are predominately affected by a growth of amplitude for decreasing values of the damping and increasing prescribed frequency. The effects of damping and prescribed frequency are perceived to behave as a more normal dynamical system. With increasing values of $\xi$, the relative displacement amplitude decreases, while increasing $f_p$ increases the relative displacement amplitude.

Table 15 details the maximum and minimum values for motion of cases A and D under the $R_{opt}\approx 1/C_p{\omega }_p$ condition. Case A motion has been shown to be extremely impactful on the motion and harvested power. The $K$ term causes the relative motion to be larger when considered to be the defined $2\pi /L$ than when zero. This undulatory term causes the harvested power to be exceptionally high with values varying with an exponential power difference compared to when it is considered zero. The relative amplitude of the head patch is smaller for case A and larger for case D motion when compared to the tail patch values. The tail patch harvests more power than the head patch for cases A and D.

The comparison of results for the uncoupled versus coupled BCFEH show that K has a large impact on the system. These systems harvest the largest amount of energy when $f_p=8 \mathrm{Hz}$. For case A, the head patch harvests $6.5542\times {10}^{-2} \mathrm{W}$ and $2.6159\times {10}^{-3} \mathrm{W}$ for the coupled and uncoupled systems, and the coupled and uncoupled tail patch harvests $6.9594\times {10}^{-2} \mathrm{W}$ and $1.0817\times {10}^{-1} \mathrm{W}$. In the coupled system, the head and tail patch harvest similar amounts of power, but the uncoupled system has a large variability of results. The maximum power for the head patch coupled system is larger than the uncoupled case. The tail patch uncoupled power is larger than the coupled system. For case D, the coupled and uncoupled head patch harvests $6.0246\times {10}^{-3} \mathrm{W}$ and $1.4288\times {10}^{-3} \mathrm{W}$, and the tail patch harvests $6.0971\times {10}^{-3} \mathrm{W}$ and $1.4288\times {10}^{-3} \mathrm{W}$. The coupled power is greater than the uncoupled system for both patches undergoing case D motion.

6. Conclusions

The development of a coupled electromechanical model for a body caudal fin energy harvester undergoing a forced body undulatory–oscillation actuation presents challenges that require multiple investigations. The undulatory controlling term influences the modal summation, in combination with terms that dramatically influence and affect the coupled system. The results show that many modes are required to get a semblance of convergence. Seven modes were used in this study, as mathematical challenges appeared when determining the mode shape solution for higher mode numbers. The Gauss law component of the coupled system has a forcing term depending on the motion, which is atypical for conventional energy harvesting systems. For high values of the structural damping term, the Gauss law component can be rearranged into a previously used uncoupled model that is used to investigate only the potential voltage harvesting of the motion.

Comparing the results of the coupled and uncoupled representations in terms of their energy harvesting potential for different cases of motion gave us an insight into how these systems could react to different configurations. The importance of this investigation lies in the identification that a structure that is successfully designed to move with a forced undulatory motion may be susceptible to additional effects resulting from the motion interacting with the structure. By removing the impactful undulating term, the performances of the patch at the head and tail locations become more alike as the complexity is reduced and the motion at these locations becomes more similar. The adaptation of a high damping architectonic approach could decouple the system and lead to a pure forced motion–electrical response that is more predictable; this could be recommended for pursuit as a possible method to reduce complexities in these prescribed actuation systems. The coupled system considering low damping leads to the coupling of the prescribed motion and structural response, which is very complex even when out of resonance. This creates a limitation of scope. If the system were designed in a way that the natural harmonics of the structure would become more influential, it would be possible for the motion–structure coupling to be damaging.