Energy harvesting using an inverted piezohydroelastic flag with resistor–inductor–capacitor circuit

Abstract

We develop a high-fidelity multi-physics model and, by using it, study the energy conversion process of a piezohydroelastic flag. Instead of a regular flag with a clamped upstream leading edge, we use an inverted flag so as to make the best of fluid elastic instability for energy harvesting. Moreover, different from many previous studies where a resistor–capacitor (RC) circuit is usually used, we adopt a resistor–inductor–capacitor (RLC) circuit for electricity generation. The influences of several key parameters associated with fluid, structure and electric dynamics are studied. Significantly different response modes are identified, among which the symmetric- and asymmetric-flutter modes are most suitable for sustainable energy harvesting, both emerging with moderate bending stiffness. If only deploying the RC circuit, increasing the resistance makes the flag more stable. By adding an inductor to turn the RC circuit into an RLC one, we observe the occurrence of ‘lock-in’ between the flag frequency and the circuit frequency as first reported on a regular flag by Xia et al. (Phys. Rev. Appl., vol. 3, 2015, 014009). This phenomenon can significantly enhance the energy output, but it only happens when the circuit resistance is sufficiently large. From a derivation based on dynamic mode decomposition analysis, we further identify an optimal condition for maximizing the energy output, which can serve as a guideline to determine whether deploying an inductor can boost the performance, and, if yes, the required inductance. The findings from this study can better guide the design of flow-induced-vibration-based piezoelectric energy harvesters for microelectronic devices.

1. Introduction

Low-power electronics have been widely used in many areas such as implanted biomedical devices (Mitcheson et al. 2008; Paulo & Gaspar 2010) and wireless sensory nodes in remote or hostile locations for environmental monitoring (Shaikh & Zeadally 2016; Kanoun 2018). Conventionally, many of these electronics are powered by chemical batteries which necessitate periodic replacement and costly maintenance. To circumvent this situation, sustainable power generation can be achieved by extracting ambient energy. Some possible ambient energy sources are, for instance, thermal energy, light energy or flow energy. The present research focuses on the extraction of a ubiquitous source of flow energy via flow-induced vibration (FIV). Once triggered by fluid elastic instability, FIV is able to induce unsteady loads on structures. Classical FIV models include fluttering airfoils (Zhu & Peng 2009; Zhu 2011; Bibo & Daqaq 2013; Wang et al. 2020; Tamimi et al. 2022) and thin plates (Sodano, Inman & Park 2004; Dunnmon et al. 2011; Giacomello & Porfiri 2011; Xu-Xu, Barrero-Gil & Velazquez 2016; Qadri, Zhao & Tang 2020; Zhao et al. 2021) mounted in a uniform flow, cantilever beams placed in the Kármán vortex street behind a cylinder (Adhikari, Rastogi & Bhattacharya 2020; Du et al. 2022), galloping of bluff bodies (Barrero-Gil, Alonso & Sanz-Andres 2010; Tan, Zuo & Yan 2021), etc.

To convert the mechanical energy from ambient flows into electricity, there are several transduction mechanisms such as piezoelectric, electromagnetic and electrostatic generators, among which the piezoelectric transducer is of the major interest over the past decade due to its inherent electromechanical coupling and transduction capacity, ease of applications and high power density (Stojcev, Kosanovic & Golubovic 2009). Piezoelectricity refers to the phenomenon that piezoelectric material becomes polarized due to the accumulation of electric charge on the electrodes in response to applied mechanical stress (Cross 2004; Panda 2009). Through this mechanism, FIV-based piezoelectric transducers are able to provide stable power. As such, they have been considered as a promising alternative to batteries for low-power electronics (Uchino 2018; Hamlehdar, Kasaeian & Safaei 2019; Safaei, Sodano & Anton 2019).

Many experimental, theoretical and numerical efforts have been spent to explore the energy conversion process of FIV-based piezoelectric transducers, on which Hamlehdar et al. (2019) have done a comprehensive review. Experimental studies (e.g. Dunnmon et al. 2011; Hobbs & Hu 2012; Song et al. 2015; Petrini & Gkoumas 2018) have demonstrated the feasibility and reliability of the FIV-based piezoelectric transducer under various conditions. Theoretical studies, such as those conducted by Abdelkefi, Nayfeh & Hajj (2012) and Tan, Yan & Hajj (2016), used decoupled electromechanical models to analyse the dynamics of FIV-based piezoaeroelastic systems. However, these decoupled models oversimplified the sophisticated fluid–structure interactions. To address this issue, a number of coupled fluid–structure–electric models have been developed in the past decade. For example, Doaré & Michelin (2011) and Michelin & Doaré (2013) developed such a multi-physics model to study a two-dimensional elastic plate placed in an axial flow. The plate experiences self-sustained flutter when the flow velocity exceeds a critical value (Connell & Yue 2007; Alben & Shelley 2008). Electricity is then generated via distributed piezoelectric patches that are applied along the plate's surface. They found that energy harvesting efficiency of this system is optimized when the flutter frequency is tuned to the time scale of the electric system. De Marqui, Erturk & Inman (2010) presented a piezoaeroelastic model combining an unsteady vortex-lattice model and an electromechanically coupled finite-element model to investigate the electrical power output of a cantilever plate placed on the flexible wing of micro air vehicles.

Although revealing useful physical insights in the energy conversion mechanism, the above coupled fluid–structure–electric models employed either the slender body theory or the potential flow theory to describe the fluid dynamics, which may not be applicable to low- or intermediate-Reynolds-number flows in which most low-power electromechanical devices are operated (Taylor et al. 2001; Wang & Ko 2010; Wang & Liu 2011). To address this issue, Akcabay & Young (2012) and Shoele & Mittal (2016) developed new sets of fully coupled fluid–structure–electric models, in which the fluid dynamics is described by the Navier–Stokes equations and the interactions between the viscous flow and the piezoelectric structure are dealt with by the immersed boundary method (Peskin 2002). In these studies, the impacts of many parameters such as density ratio, Reynolds number, electromechanical coupling coefficient, etc. upon the piezohydroelastic system have been extensively investigated. The difference of these two studies lies in the piezoelectric patch configuration: the former utilized the single-patch configuration while the latter used distributed piezoelectric patches.

Existing studies have also revealed some scenarios where the energy harvesting performance of FIV-based piezoelectric transducers can be enhanced. Such improvements can be achieved through manipulations in fluid, structure or electrical domains, either separately or collectively. For example, from the fluid perspective, Mazharmanesh et al. (2022) found that, compared with in uniform flows, a piezoelectric flag placed in oscillating flows is able to create much higher energy due to constructive interactions between the vorticity flow and the flag. From the structure perspective, inverted flags, clamped at the downstream trailing edge and free to flap at the upstream leading edge, are found to be beneficial to energy harvesting (Kim et al. 2013; Ryu et al. 2015; Tang, Liu & Lu 2015). Compared with regular flags with clamped upstream leading edge, inverted flags are far more unstable and are able to produce sustained large-amplitude vibrations in a wider range of bending rigidity and flow speeds (Shoele & Mittal 2016).

As for the electrical circuit, most existing studies only considered the simplest resistor–capacitor (RC) circuit, i.e. a capacitor rendered by the piezoelectric material's intrinsic feature and a resistor representing the electrical load (Doaré & Michelin 2011; Michelin & Doaré 2013; Dias, De Marqui & Erturk 2014; Shoele & Mittal 2016; Safaei et al. 2019; Mazharmanesh et al. 2022). Energy harvesting performance can be enhanced by introducing an inductor into this conventional RC circuit, forming a resistor–inductor–capacitor (RLC) circuit that intrinsically has a resonant frequency. By inheriting the fully coupled model used in Michelin & Doaré (2013) and deploying an RLC circuit on a regular flag, Xia, Michelin & Doaré (2015) found that the flapping frequency of the flag was locked onto the circuit's natural frequency under certain circumstances, leading to a wider operating frequency bandwidth and an enhanced efficiency. A similar observation was also made by Wang et al. (2016). Moreover, by applying a decoupled model to a regular flag, Tan & Yan (2017) concluded that simultaneous tuning of the resistance and the inductance can achieve better energy harvesting efficiencies than tuning the resistance alone.

Clearly, the existing studies have suggested that inverted flags can generally perform better than regular flags in terms of energy generation, and RLC-circuit-based flags may perform better than RC-circuit-based flags. Intuitively, one would expect that RLC-circuit-based inverted flags can exhibit the best performance. However, so far, no study has been reported on such a highly promising configuration. To this end, the current study aims at developing a fully coupled fluid–structure–electric model and using it to study the energy conversion of an inverted piezohydroelastic flag mounted in a viscous flow. The fluid dynamics will be depicted with the Navier–Stokes equations and the fluid–structure coupling will be handled by the immersed boundary method. An RLC circuit will be applied for the electricity generation. With this multi-physics model, we explore the characteristics of the inverted piezohydroelastic flag system, particularly the influence of the RLC circuit on the dynamics of the flag and the resulting energy harvesting performance.

2. Physical problem

As shown in figure 1(a), the inverted flag we consider here is reduced to a cantilevered plate with length $L$ and infinite span, subjected to a parallel fluid flow of speed $U_\infty$ moving from the free leading edge to the clamped trailing edge. The plate is assumed to be thin and inextensible, and only its out-of-plane deformation is considered. As such, we can use a two-dimensional Euler–Bernoulli beam model to describe its structural dynamics.

Figure 1. (a) Schematic of the inverted piezoelectric flag. (b) Electric circuit of a piezoelectrical patch pair connected with the output loop and the equivalent RLC circuit.

The lower and upper surfaces of the flag are covered by piezoelectric layers (referred to as the bimorph configuration (Bonello & Rafique 2011); see figure 1b), which contain continuous piezoelectrical patch pairs of infinitesimal length (Bisegna, Caruso & Maceri 2006; Doaré & Michelin 2011). The patch pairs are connected in series with opposite polarity. The electrodes are coupled with the power output circuit consisting of a resistor and an inductor, which are connected in parallel. As the piezoelectric material is deformed by the local curvature of the flag, the electrodes will have a potential difference so that an electric current can be generated in the circuit. The power dissipated by the resistor as the electric current passes through is the harvested energy.

3. Problem formulations

3.1. Governing equations for the coupled fluid–structure–electric system

Described in the Eulerian coordinate $\boldsymbol {x}\equiv (x,y)$, the fluid dynamics is governed by the incompressible Navier–Stokes (N–S) equations with constant fluid density $\rho$ and viscosity $\mu$. The dimensionless N–S equations are

(3.1)\begin{equation} \left. \begin{gathered} \frac{\partial \boldsymbol{u}}{\partial t}+\boldsymbol{\nabla}\boldsymbol{\cdot}(\boldsymbol{u}\boldsymbol{u})=-\boldsymbol{\nabla} p+\frac{1}{Re}\nabla^2\boldsymbol{u}+\boldsymbol{f}(\boldsymbol{x})\quad \text{and}\\ \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{u}=0, \end{gathered} \right\} \end{equation}

in which $\boldsymbol {u}$ represents the fluid velocity vector, $p$ the pressure, $Re\equiv \rho U_\infty L/\mu$ the Reynolds number and $\boldsymbol {f}$ the fluid–structure interaction force density.

The three-layer ‘sandwich’ in figure 1 is modelled as an Euler–Bernoulli beam. Let $\boldsymbol {X}(s,t)$ denote the instantaneous displacement vector of an arbitrary structural element at the Lagrangian location $s$, the conservation of momentum and the inextensibility constraint (i.e. $\partial \boldsymbol {X}/\partial s\boldsymbol {\cdot }\partial \boldsymbol {X}/\partial s=1$) for the flag lead to

(3.2) \begin{equation} \left. \begin{gathered} m_s\frac{\partial^2\boldsymbol{X}(s,t)}{\partial t^2} =\frac{\partial}{\partial s}\left[\sigma(s,t)\boldsymbol{e}_\tau-\frac{\partial M(s,t)}{\partial s}\boldsymbol{e}_n\right]-\boldsymbol{F}(s,t)+ m_s\boldsymbol{g},\\ \boldsymbol{X} =\boldsymbol{X}_0 \quad \text{and} \quad \frac{\partial\boldsymbol{X}}{\partial s}=(1,0), \quad \text{at the trailing edge}\\ M =0 \quad \text{and} \quad \sigma\boldsymbol{e}_\tau-\frac{\partial M}{\partial s}\boldsymbol{e}_n=0, \quad \text{at the leading edge} \end{gathered} \right\} \end{equation}

where $m_s$ is the mass per unit length of the beam and $\boldsymbol {F}$ stands for the hydrodynamic loading. The gravity acceleration $\boldsymbol {g}$ is set as zero throughout the manuscript unless otherwise stated (it is only involved in Appendix A.1). Here, $\boldsymbol {e}_n$ and $\boldsymbol {e}_\tau$ are, respectively, the unit normal and tangential vectors, defined as $\boldsymbol {e}_\tau =\partial \boldsymbol {X}/\partial s$ and $\boldsymbol {e}_n=(0,0,1)\times \boldsymbol {e}_\tau$ and $\boldsymbol {X}_0$ is the initial position of the flag. The tension force $\sigma$ is determined by the constraint of inextensibility(Huang, Shin & Sung 2007)

(3.3)\begin{align} \frac{\partial^2}{\partial s^2}(\sigma\boldsymbol{e}_\tau)\boldsymbol{\cdot}\boldsymbol{e}_\tau=\frac{m_s}{2}\frac{\partial^2}{\partial t^2}(\boldsymbol{e}_\tau\boldsymbol{\cdot}\boldsymbol{e}_\tau)-m_s\frac{\partial\boldsymbol{e}_\tau}{\partial t}\boldsymbol{\cdot}\frac{\partial\boldsymbol{e}_\tau}{\partial t}+\frac{\partial}{\partial s}\left[\frac{\partial}{\partial s}\left(\frac{\partial M}{\partial s}\boldsymbol{e}_n\right)+\boldsymbol{F}-m_s\boldsymbol{g}\right]\boldsymbol{\cdot}\boldsymbol{e}_\tau . \end{align}

In (3.2), $M$ is the total internal bending moment of the piezoelectric flag resulting from the flag's deformation and the reverse piezoelectric effect of the piezoelectric patches, i.e.

(3.4)\begin{equation} M=k_b\frac{\partial^2\boldsymbol{X}}{\partial s^2}\boldsymbol{\cdot}\boldsymbol{e}_n-\chi\upsilon(s,t), \end{equation}

in which $k_b$ is the bending stiffness, $\chi$ the structure–electric coupling coefficient and $\upsilon$ the electric voltage. We consider a continuous distribution of patches (Doaré & Michelin 2011; Michelin & Doaré 2013) such that the voltage $\upsilon$ and electrical charge displacement $Q$ are continuous functions of the flag's location $s$. Based on Gauss’ law (Jackson & David 1998), the electrical equation for the RLC circuit is

(3.5)\begin{equation} \left. \begin{gathered} \frac{\partial Q}{\partial t}+\gamma\upsilon+\int_0^t\frac{\upsilon(\tau)}{\ell}\,{\rm d}\tau=0,\\ Q=c\upsilon+\chi\frac{\partial^2\boldsymbol{X}}{\partial s^2}\boldsymbol{\cdot}\boldsymbol{e}_n. \end{gathered} \right\} \end{equation}

It can be further reduced to

(3.6)\begin{equation} c\frac{\partial^2 \upsilon}{\partial t^2}+\gamma\frac{\partial \upsilon}{\partial t}+\frac{\upsilon}{\ell}=-\chi\frac{\partial^2}{\partial t^2}\left(\frac{\partial^2\boldsymbol{X}}{\partial s^2}\boldsymbol{\cdot}\boldsymbol{e}_n\right) , \end{equation}

where $c$ is the intrinsic capacitance of the piezoelectric patch pair, $\gamma$ the lineic conductivity of the resistor (inverse of the resistance $R$) and $\ell$ the inductance.

The coupling of the fluid and piezoelastic dynamics is accomplished with the immersed boundary method (IBM). The IBM-based fluid–structure interaction algorithm has been well documented in past studies (Peskin 2002; Bi & Zhu 2019), herein, we only briefly sketch its key points. The fluid–structure interacting forces $\boldsymbol {f}(\boldsymbol {x},t)$ and $\boldsymbol {F}(s,t)$ are defined in the Eulerian and Lagrangian coordinates, respectively. The transformation between these two quantities is achieved by using the Dirac delta function

(3.7)\begin{equation} \boldsymbol{f}(\boldsymbol{x},t) =\int_\varGamma\boldsymbol{F}(s,t)\delta(\boldsymbol{x}-\boldsymbol{X}(s,t))\,{\rm d}s, \end{equation}

in which $\varGamma$ represents the Lagrangian domain. Following Goldstein, Handler & Sirovich (1993), $\boldsymbol {F}$ is calculated by

(3.8)\begin{equation} \boldsymbol{F}(s,t)=\tilde{\alpha}\int_0^t\left[\boldsymbol{U}(s,\tau) -\boldsymbol{V}(s,\tau)\right]\,{\rm d}\tau+\tilde{\beta}\left[\boldsymbol{U}(s,t)-\boldsymbol{V}(s,t)\right], \end{equation}

where $\tilde {\alpha }$ and $\tilde {\beta }$ are sufficiently large negative constants (Goldstein et al. 1993), $\boldsymbol {V}(s,t)$ is the structural velocity, i.e. $\boldsymbol {V}=\partial \boldsymbol {X}/\partial t$, $\boldsymbol {U}(s,t)$ denotes the local fluid velocity at the immersed structural position $s$, which is calculated by

(3.9)\begin{equation} \boldsymbol{U}(s,t)=\int_\varOmega\boldsymbol{u}(\boldsymbol{x},t)\delta(\boldsymbol{x}-\boldsymbol{X}(s,t))\,{\rm d}\kern0.7pt \boldsymbol{x}, \end{equation}

where $\varOmega$ denotes the fluid domain. Since the fluid and structure dynamics are solved independently, the matching of $\boldsymbol {U}$ and $\boldsymbol {V}$ is achieved by (3.8), through which a penalty force will be generated, which converges to the physical fluid–structure interaction force after iterations.

The energy transfer pathway of the piezohydroelastic system is presented in figure 2. The ambient flow powers the flag at a rate of $P_{in}$, leading to substantial mechanical energy $\mathscr {E}_m$ (including the kinetic energy and the strain energy) being stored in the flag. A portion of $\mathscr {E}_m$ is converted to electric energy $\mathscr {E}_{e}$ at a rate of $P_{me}$, which is stored in the piezoelectric capacitor and the inductor in the output electric circuit. The electricity will then be dissipated by the resistor at a rate of $P_r$; $P_{r}$ is always positive according to its nature, but the instantaneous values of the other two energy transfer rates $P_{in}$ and $P_{me}$ can be negative during certain time intervals. For this reason, reversed arrows are added beside the main arrows. Nevertheless, the energy is mainly transferred along the main arrows and the time-averaged values of these two energy transfer rates are positive. According to the law of energy conservation, we have

(3.10a,b)\begin{equation} \frac{{\rm d}\mathscr{E}_m}{{\rm d}t}=P_{in}-P_{me}, \quad \frac{{\rm d}\mathscr{E}_{e}}{{\rm d}t}=P_{me}-P_{r}. \end{equation}

When the system achieves its periodic steady state, we have $\langle P_r\rangle =\langle P_{me}\rangle =\langle P_{in}\rangle$, where $\langle {\cdot }\rangle$ denotes the time-averaging operator.

Figure 2. The energy transfer pathway of a piezohydroelastic flag.

The performance of the system is evaluated through energy harvesting efficiency $C_p$, the ratio between the time-averaged power expenditure by the resistor and the theoretical upper limit of the flow-energy flux past the projected area that is swept by the fluttering flag, i.e.

(3.11)\begin{equation} C_p=\frac{\langle P_r\rangle}{\dfrac{1}{2}\rho 2L U_\infty U_\infty^2}=\frac{\langle P_r\rangle}{\rho LU_\infty^3}. \end{equation}

In the present work, the characteristic length, velocity, density and voltage are chosen as $L$, $U_\infty$, $\rho$ and $U_\infty \sqrt {\rho L/c}$, respectively. The electroelastic equations ((3.2) and (3.6)) can be non-dimensionalized as

(3.12)\begin{equation} \left. \begin{gathered} \frac{1}{m^*}\frac{\partial^2\boldsymbol{X}}{\partial t^2 }=\frac{\partial}{\partial s}\left(\sigma(s,t)\boldsymbol{e}_\tau\right)- \frac{\partial}{\partial s}\left[K^*\frac{\partial^3\boldsymbol{X}}{\partial s^3}-\alpha\sqrt{K^*}\frac{\partial \upsilon}{\partial s}\boldsymbol{e}_n\right]-\boldsymbol{F}(s,t),\\ \frac{\partial^2\upsilon}{\partial t^2} +\frac{1}{\beta}\frac{\partial\upsilon}{\partial t}+\omega_0^2\upsilon=-\alpha \sqrt{K^*}\frac{\partial^2}{\partial t^2}\left(\frac{\partial^2\boldsymbol{X}}{\partial s^2}\boldsymbol{\cdot}\boldsymbol{e}_n\right), \end{gathered} \right\} \end{equation}

in which $m^*$ represents the inertial ratio between the fluid and structure, $K^*$ the normalized bending stiffness, $\alpha$ the normalized electromechanical coupling coefficient, $\beta$ the normalized resistance and $\omega _0$ the natural frequency of the electrical circuit, defined as follows:

(3.13a–e)\begin{equation} m^*=\frac{\rho L}{m_s}, \quad K^*=\frac{k_b}{\rho {U_\infty}^2 L^3}, \quad \alpha=\frac{\chi}{\sqrt{k_bc}}, \quad \beta=\frac{cU_\infty}{\gamma L}, \quad \omega_0=\frac{L}{U_\infty\sqrt{\ell c}}. \end{equation}

Note that Michelin & Doaré (2013) and Shoele & Mittal (2016) used $1/{U^*}^2$ to replace $K^*$, where $U^*$ is defined as the normalized free-stream velocity. In addition, the power harvesting efficiency $C_p$ can also be viewed as the non-dimensional time-averaged energy harvested by the system:

(3.14)\begin{equation} C_p=\left.\left\langle\int_0^1\upsilon^2\,{\rm d}s\right\rangle\right/\beta. \end{equation}

Hereafter, all simulation results will be presented in the non-dimensional form.

3.2. Numerical implementation

The hydrodynamic equations (3.1) are spatially discretized using the finite difference approach over a set of staggered Cartesian grids. The Crank–Nicolson scheme is used to temporally discretize the diffusion and convection terms. The discrete form of the equations becomes

(3.15)\begin{align} \left. \begin{gathered} \frac{\boldsymbol{u}^{n+1}-\boldsymbol{u}^n}{\Delta t}+\frac{\mathscr{H}(\boldsymbol{u}^{n+1})+\mathscr{H}(\boldsymbol{u}^n)}{2}=-\mathscr{G}(\,p^{n+1/2})+\frac{\mathscr{L}(\boldsymbol{u}^{n+1})+\mathscr{L}(\boldsymbol{u}^n)}{2Re}+\boldsymbol{f}^{n+1},\\ \mathscr{D}(\boldsymbol{u}^{n+1})=0, \end{gathered} \right\} \end{align}

where $\Delta t$ is the time step, the superscript ‘$n+1$’ denotes the $(n+1)$th time step. Here, $\mathscr {H}$ represents the discrete convective operator, $\mathscr {L}$ the discrete Laplace operator, $\mathscr {G}$ the discrete gradient operator and $\mathscr {D}$ the discrete divergence operator. Following Kim, Baek & Sung (2002), the discrete hydrodynamic equations are solved with a projection method. The general process of updating the flow velocity and pressure fields is well documented in this reference.

Regarding the numerical treatment for the piezoelastic equations, uniform staggered grids are also used to discretize the body along $s$. A fully implicit time advancement method is employed, leading to the discrete piezoelastic equations

(3.16) $$\begin{align} \frac{1}{m^*}\frac{\boldsymbol{X}^{n+1}-2\boldsymbol{X}^{n}+\boldsymbol{X}^{n-1}}{\Delta t^2}&=D_s(\sigma^{n+1/2}D_s\boldsymbol{X}^{n+1})-K^*D_{ssss}\boldsymbol{X}^{n+1}\nonumber\\ &\quad +\alpha\sqrt{K^*}D_s(D_s\upsilon^{n+1}\boldsymbol{e}_n^{n+1})-\boldsymbol{F}^{n+1}, \end{align}$$

and

(3.17) $$\begin{align} &\frac{\upsilon^{n+1}-2\upsilon^n+\upsilon^{n-1}}{\Delta t^2}+\frac{\upsilon^{n+1}-\upsilon^n}{\beta\Delta t}+\omega_0^2\upsilon^{n+1}\nonumber\\ &\quad =\alpha\sqrt{K^*}\frac{D_{ss}\boldsymbol{X}^{n+1}\boldsymbol{\cdot}\boldsymbol{e}_n^{n+1}-2D_{ss}\boldsymbol{X}^{n}\boldsymbol{\cdot}\boldsymbol{e}_n^n+D_{ss}\boldsymbol{X}^{n-1}\boldsymbol{\cdot}\boldsymbol{e}_n^{n-1}}{-\Delta t^2}, \end{align}$$

where $D_s$, $D_{ss}$ and $D_{ssss}$ are the first-, second- and fourth-order central difference operators with respect to $s$, respectively. The tension $\gamma ^{n+1/2}$ is calculated by

(3.18)$$\begin{align} D_{ss}(\sigma^{n+1/2}\boldsymbol{e}_\tau^\prime)\boldsymbol{\cdot}\boldsymbol{e}_\tau^\prime&=\frac{1-2(\boldsymbol{e}_\tau\boldsymbol{\cdot}\boldsymbol{e}_\tau)^n+(\boldsymbol{e}_\tau\boldsymbol{\cdot}\boldsymbol{e}_\tau)^{n-1}}{2m^*\Delta t^2}-\frac{1}{m^*}(D_s\boldsymbol{V}^\prime\boldsymbol{\cdot} D_s\boldsymbol{V}^\prime)\nonumber\\ &\quad +D_s\left[K^*D_{ssss}\boldsymbol{X}^\prime-\alpha\sqrt{K^*}D_s(D_s\upsilon^\prime\boldsymbol{e}_n^n)+\boldsymbol{F}^{\prime}\right]\boldsymbol{\cdot}\boldsymbol{e}_\tau^\prime, \end{align}$$

in which the superscript prime represents the intermediate time step, indicating the variables will be updated during feedback iterations.

The following time-marching method including two feedback loops is chosen to resolve the coupled hydrodynamic and piezoelastic equations:

1. (i) Solve (3.15) for $\boldsymbol {u}^{n+1}$.

2. (ii) Update the penalty force $\boldsymbol {F}$ via (3.8) and (3.9) by using the updated $\boldsymbol {u}^{n+1}$.

3. (iii) Solve (3.16) and (3.18) for $\boldsymbol {X}^{n+1}$, then update $\boldsymbol {F}$ through (3.8)–(3.9).

4. (iv) Update $\upsilon ^{n+1}$ via (3.17).

5. (v) Repeat steps (iii) to (iv) until the convergence of the piezoelastic solutions.

6. (vi) Update $\boldsymbol {f}$ by using (3.7).

7. (vii) Repeat steps (i) to (vi) until the fluid field $\boldsymbol {u}^{n+1}$ converges.

More details of this numerical scheme can be found in Shoele & Mittal (2016).

All simulations will be performed in the rectangular domain depicted in figure 3, where the boundary conditions are also specified. There are $500\times 400$ grid nodes in the $x$ and $y$ directions. A uniform grid size $\Delta x=\Delta y=0.011$ is used in the vicinity of the flag, whose structural grid size is chosen to be $\Delta s=0.0067$. A constant time step $\Delta t=0.0002$ is employed. All these numerical parameters are determined based on a series of sensitivity tests.

Figure 3. Computational domain (not to scale) and boundary conditions.

The accuracy of our fluid–structure–electrical numerical model is confirmed multi-hierarchically through comparisons with benchmark results (see Appendix A).

4. Results and discussion

We conduct multi-physics simulations using the above coupled structure–fluid–electrical model to study the energy harvesting performance of the inverted piezohydroelastic flag. The dynamics of the system is dependent upon the five non-dimensional parameters listed in (3.13a–e) plus the Reynolds number $Re$. As revealed in Shoele & Mittal (2016), the dynamics of the piezoelectric flag is very similar when $Re$ lies in the range of $25\unicode{x2013} 800$. Hence, in the present study $Re$ is fixed at 200, which is sufficiently high to capture the key features of the system and, in the meantime, not computationally expensive. Although the density ratio, $m^*$, is also a key factor determining the dynamics and energy harvesting performance of the system, its influence has been extensively studied in literature and hence is fixed at $m^*=1.0$ in the present study. The electromechanical coupling coefficient, $\alpha$, describes the coupling strength between the mechanical part and the output circuit. In this study, we only focus on the strong coupling scenario in favour of high energy conversion rate, and choose a high value, i.e. $\alpha =0.9$, unless otherwise mentioned. As such, we will concentrate on the influences of the remaining three parameters, i.e. $K^*$, $\beta$ and $\omega _0$, respectively characterizing the structural bending stiffness, resistive and inductive properties of the electrical circuit, on the dynamics and energy harvesting performance of the system.

4.1. Structural response

In what follows, we explore the structural dynamics of the inverted piezohydroelastic flag with three representative bending rigidity values $K^*=0.1$, 0.25, 0.4 (0.1 represents a super flexible flag while 0.4 connotes a very stiff one), in a large electrical parameter space, i.e. $10^{-2}\le \beta \le 10^2$ and $0\le \omega _0\le 10$. Since our simulations show that the flag tends to stay undeformed (referred to as the ‘straight mode’) (Kim et al. 2013) as the bending rigidity reaches the level of $K^*=0.5$, this scenario will be excluded from our consideration.

Four dynamic modes are observed, namely the symmetric-flutter mode where the flag performs symmetric-flutter motion about its equilibrium position, the small-deflection mode where the flag slightly deforms to one side, the large-deflection mode where the flag significantly deforms to one side and the asymmetric-flutter mode where the flag performs flutter motion but the upward and downward strokes are asymmetric about its equilibrium position. Figure 4 presents the time histories of the lateral displacement of the flag free tip, $y_{tip}$, and the body profiles of these four dynamic modes. These modes are consistent with those observed in previous studies (Kim et al. 2013; Cerdeira et al. 2021), in which the dynamics of inverted flags (without electric circuits) was experimentally investigated for various bending rigidity values and mass ratios, indicating that the inclusion of electric circuits does not induce any new response mode. Only when the two flutter modes are exhibited, i.e. the symmetric- and asymmetric-flutter modes, does the flag undergo large-amplitude oscillations so that it can continuously and effectively extract flow energy.

Figure 4. (a–d) Deformation and wake patterns of various dynamic modes: (a) the symmetric-flutter (SF) mode, (b) the small-deflection (SD) mode, (c) the large-deflection (LD) mode, (d) the asymmetric-flutter (AF) mode. The contours are based on vorticity (dashed contours represent counterclockwise vortices with negative vorticity). The insets present the corresponding full body profiles. (e) Time histories of the vertical displacement of the flag tip in different dynamic modes. Panels show (a) $K^* = 0.25$, $\beta = 0.1$, $\omega _0 = 10$; (b) $K^* = 0.4$, $\beta = 1.0$, $\omega _0 = 0.5$; (c) $K^* = 0.1$, $\beta = 10$, $\omega _0 = 5.0$; (d) $K^* = 0.25$, $\beta = 10$, $\omega _0 = 0$.

Figure 4 also presents the wake patterns of these dynamic modes. As shown in figure 4(b), the wake of the small-deflection mode is characterized by a pair of shear layers, developed from both the upper and lower surfaces of the flag, which are elongated and eventually break up due to the Kelvin–Helmholtz instability. In the large-deflection mode, the large mean deformation and small-amplitude oscillation make the flag behave like a bluff body. As a result, vortices of opposite signs alternatively shed from the leading and trailing edges of the flag, leading to a typical Kármán vortex street in the wake (i.e. the S mode) as shown in figure 4(c). In the symmetric-flutter mode shown in figure 4(a), a pair of opposite-signed vortices shed from the leading edge once the flag reaches either one of its stroke extremes, which then propagate downstream with a prominent angle. As a result, two streets of vortex pairs (i.e. the 2P mode) appear in the wake. Differently, in the asymmetric-flutter mode the flag flutters with a much smaller amplitude and a lower frequency, as can be read in figure 4(e). As such, only one pair of vortices are generated during one flapping period, forming a 2S-mode vortex street with very small lateral distance, as shown in figure 4(d).

Figure 5 maps the structural dynamics of the inverted piezohydroelastic flag in the electrical parameter space. Here, we use green markers (both up-pointing and down-pointing triangles) to denote the symmetric- and asymmetric-flutter modes that are appropriate for flow-energy harvesting. When the flag is super flexible with $K^*=0.1$, there is only a small green region in the upper left corner of figure 5(a), while all the remaining cases belong to the unsuitable large-deflection mode. This is not surprising because a too soft flag fails to provide sufficient recovery bending moment for sustainable flutter motion. On the other hand, as the flag becomes very stiff with $K^*=0.4$ (figure 5c), the cases in the upper left corner are in the small-deflection mode or even remain undeformed (i.e. in the straight mode), and hence inappropriate for flow-energy harvesting. Only when the flag is moderately flexible with $K^*=0.25$ does the system undergo the symmetric- or asymmetric-flutter motions in most of the map except for the case with the small-deflection mode at $\beta =10^2$ and $\omega _0=0$ (see figure 5b).

Figure 5. Distribution of dynamic modes in the $\beta -\omega _0$ space for flags of different bending stiffness. Panels show (a) $K^* = 0.1$, (b) $K^* = 0.25$, (c) $K^* = 0.4$.

Note that all the cases in the first column of each panel of figure 5 are only equipped with an RC circuit, where $\omega _0=0$. Furthermore, with the smallest resistance (i.e. $\beta =10^{-2}$) the bottom left case approximates the short-circuit case (or, equivalently, the case without piezoelectric effect, i.e. $\alpha =0$), whereas with the largest resistance (i.e. $\beta =10^2$), the upper left case approximates the open-circuit case. From the results in this column, it is interesting to see that increasing the resistance in the output RC circuit makes the flag appear stiffer, or in other words, more stable. Specifically, as $\beta$ increases, the flag changes from the large-deflection mode to the symmetric-flutter mode in figure 5(a), from the symmetric-flutter mode to the asymmetric-flutter mode and then to the small-deflection mode in figure 5(b) and from the symmetric-flutter mode all the way to the straight mode in figure 5(c).

Figure 6 further shows the variations of the flag's dominant flapping frequency $\omega$ and amplitude $A$ against the electrical circuit's undamped natural frequency $\omega _0$ at selected resistance and bending stiffness values. Here, $\omega$ is calculated by Fourier analysis of the time histories of vertical displacement of the flag tip, $y_{tip}$, and $A$ is obtained by measuring the difference between the upper and lower bounds of $y_{tip}$. Since in the present study the flag neither experiences high-order oscillations nor reverses its leading-edge tip, $A$ is positively correlated with its maximum curvature. Mode types identified in figure 5 are also marked in the figure with the same set of symbols. As baseline values, the frequency and amplitude of the flag with no electromechanical coupling (i.e. $\alpha =0$), referred to as $\omega _{s}$ and $A_{s}$, respectively, are also presented as black solid straight lines.

Figure 6. Variations of flapping frequency (a,c) and amplitude (b,d) against $\omega _0$ for the moderate flexible flag ((a,b), $K^*=0.25$) and the less flexible flag ((c,d), $K^*=0.4$) equipped with various resistances $\beta$ in the electrical circuit. Here, $\omega _s$ and $A_s$, indicated by black solid lines, are the flapping frequency and amplitude of the flag with no electromechanical coupling, respectively.

When the resistance is trivial (e.g. $\beta =0.01$), the flapping frequency and amplitude are very close to the baseline values, because under this condition the RLC circuit is almost short circuited by the resistor. As a result, the circuit has almost no influence on the flag deformation, which is reminiscent of the electroelastic coupling coefficient $\alpha$ being zero. Similarly, when $\omega _0$ is sufficiently large (equivalently when the lineic admittance $1/\ell$ of the inductor is sufficiently small according to (3.13a–e)), the RLC circuit is then short circuited by the inductor. Consequently, the structural responses with respect to both the flapping frequency and amplitude converge to the baseline values asymptotically regardless of the $\beta$ value.

Only when $\beta$ is not too small and $\omega _0$ is not too large does the electrical circuit exert a significant impact on the structural dynamics. It shows that varying $\omega _0$ may lead to mode transition. For example, in figure 6(a) the asymmetric-flutter mode ($\boldsymbol {\nabla }$, lime) appearing in the RC circuit is transitioned into the symmetric-flutter mode ($\Delta$, green) after implementing a small inductor ($\omega _0=0.1$) into the circuit. The mode transition in figure 6(c), however, takes place in the RLC circuit, where the most noticeable one is between the small-deflection mode ($\Diamond$, red) and the symmetric-flutter mode ($\Delta$, green).

These mode transitions are usually accompanied by an abrupt change in the flag's dominant flapping frequency $\omega$ and amplitude $A$. In particular, as shown in figure 6(a,c), $\omega$ plunges during the transition, followed by a continuous rise up to $\omega _s$ with a slope close to that of the line $\omega =\omega _0$, suggesting a ‘lock-in’ phenomenon. As revealed in the following sections, in this ‘lock-in’ regime, the flag will resonate with the circuit, leading to an enhanced energy output. A similar observation has been made by Wang et al. (2016) when examining a regular piezoelectric flag, who also revealed that during this ‘lock-in’ the flag is able to flutter at much lower speeds. Note that, in both studies, $\omega$ is consistently smaller than $\omega _0$, the undamped natural frequency of the circuit, in the ‘lock-in’ regime, different from what was reported in Xia et al. (2015), where the two frequencies coincide with each other. This frequency reduction reflects the damped and nonlinear nature of the current coupled fluid–structure–electrical system. This ‘lock-in’ becomes less pronounced as $\beta$ decreases towards the short-circuit scenario. As shown in figure 6(a,c), the ‘lock-in’ phenomenon even vanishes when $\beta \le 0.4$ for the moderate flexible flag (i.e. $K^*=0.25$) and $\beta \le 0.1$ for the less flexible flag (i.e. $K^*=0.4$).

4.2. Energy harvesting performance

The electrical parameters $\beta$ and $\omega _0$ characterize the resistive and inductive properties of the circuit, respectively, affecting the energy extraction performance of the inverted piezohydroelastic flag. In this section, we explore how they influence the energy output of the flag and pin down the underlying physics.

4.2.1. The RC circuit

The impact of $\beta$ on the energy output in the scenario of the RC circuit (i.e. $\omega _0=0$) is presented in figure 7. Since the stiff flag (i.e. $K^*=0.4$) displays multiple modes in this scenario, here, we only consider the moderate flexible flag (i.e. $K^*=0.25$) for simplicity. It is not surprising to see that trivial energy is harvested when $\beta$ is too small ($\beta \ll 1$) or too large ($\beta \gg 1$), corresponding to the short- and open-circuit conditions, respectively. A peak energy output is then achieved at a moderate $\beta$ value between 0.1 and 1. Similar observations have also been reported in Michelin & Doaré (2013) and Shoele & Mittal (2016). They attributed the optimal performance to the ‘resonance’ achieved by matching the electrical time scale and the mechanical time scale, i.e. $\omega \beta \sim O(1)$. In the present study, the energy peak appears at approximately $\omega \beta =0.22$, as denoted in figure 7, confirming the above statement on optimal condition but near the lower bound of $O(1)$. Our discussions in the following part and the result presented in figure 16 of Appendix A.3 further suggest that this near-the-lower-bound optimal condition can be attributed to the enhanced damping due to the much stronger electromechanical coupling $\alpha =0.9$ applied here (instead of $\alpha =0.5$ in Michelin & Doaré 2013; Shoele & Mittal 2016).

Figure 7. (a) Variation of the energy output $C_{p}$ and its dominant DMD mode $C_{p,1}$ against $\beta$ for the moderate flexible lag with $K^*=0.25$ and $\omega _0=0$. The insets show the flag deformation for the cases marked with coloured symbols. (b) Variation of the amplitude of the dominant DMD mode against $\beta$.

In what follows, we attempt to elaborate on the energy transfer mechanism through mode analysis. As described in (3.12), the circuit undergoes damped oscillations forced by an excitation term originated from the curvature of the flag, i.e. $\mathscr {\varLambda }(s,t)=({\partial ^2\boldsymbol {X}}/{\partial s^2})\boldsymbol {\cdot }\boldsymbol {e}_n$. A dynamic-mode-decomposition (DMD) analysis (Kutz et al. 2016) is then conducted on a time series of snapshots of $\mathscr {\varLambda }$ with an equal time interval $\Delta T=0.0125$. The electrical energy output from each DMD mode can be formulated as

(4.1)\begin{equation} C_{p,j}=\frac{2\alpha^2K^*|b_j|^2}{\beta+\dfrac{1}{\beta {\rm Im}\{\omega_j\}^2}}, \end{equation}

where $|b_j|$ stands for the amplitude of the $j$th DMD mode, $\textrm {Im}\{\omega _j\}$, i.e. the imaginary part of $\omega _j$ is the frequency of this mode. The derivation of (4.1) is provided in Appendix B.

As an example, in figure 8(a) we present the DMD spectrum of the case with $\beta =0.3$, where the energy output $C_p$ achieves its peak value (see figure 7a). It is seen that the dynamics of the excitation term of (3.12) can be represented by the dominant DMD mode, whose amplitude $|b_j|=12.2$. Its oscillation frequency $\textrm {Im}\{\omega _j\}=0.76$ is, not surprisingly, consistent with the flapping frequency of the flag which can be read from figure 6. Figure 8(b) plots the time evolution of $\varLambda$ at the mid-point ($s=0.5$) and its reconstruction using the first 10 DMD modes, proving the accuracy of the data decomposition process.

Figure 8. (a) The DMD spectrum of $\varLambda$ for the case with $\beta =0.3$. (b) Time history of $\varLambda$ at the flag mid-point $s=0.5$ and the reconstruction with the first 10 DMD modes. (c) Energy output $C_{p,j}$ contributed by the DMD modes.

Figure 8(c) further presents the energy contribution of each mode calculated by (4.1). As expected, the leading mode provides the greatest portion of the overall energy generation. Specifically, $C_{p,1}$ (here, the subscript ‘1’ refers to the leading DMD mode) accounts for $91\,\%$ of the total electrical energy output $C_p$. Cases in figure 7 with other $\beta$ values, without any exception, share a similar feature: the overall energy output is mostly attributed to the leading DMD mode. As a result, as displayed in figure 7(a), in which we present $C_p$ and $C_{p,1}$ together, the two quantities are very close to each other over the range of $\beta$ we consider.

From (4.1), one can see that the energy transfer is mainly determined by three variables, namely the resistance parameter, $\beta$, the magnitude, $|b_1|$, and the frequency, $\textrm {Im}\{\omega _1\}$, of the leading DMD mode of the curvature distribution. Shoele & Mittal (2016) reported that, at $\alpha =0.5$, where the electromechanical coupling is relatively weak, the flapping amplitude of the flag (which can be approximately reflected by $|b_1|$) remains nearly unchanged over $\beta$. Meanwhile, as revealed in figure 6, the flapping frequency $\omega$ (approximately represented by the leading frequency $\textrm {Im}\{\omega _1\}$) usually falls in a much narrower range (between $O(10^{-1})$ and $O(10^0)$) compared with $\beta$. Therefore, the optimal energy output $C_p$ (approximately represented by $C_{p,1}$) can be achieved roughly under the condition $\omega \beta \approx \textrm {Im}\{\omega _1\} \beta =1$, where the denominator of the right-hand side of (4.1) is minimized while the numerator is constant. This derivation, again, confirms the optimal condition proposed by Michelin & Doaré (2013) and Shoele & Mittal (2016), i.e. $\omega \beta \sim O(1)$. Nevertheless, the strong electromechanical coupling considered in the present study makes $|b_1|$ decrease prominently with $\beta$, especially in range of $10^{-1}<\beta <10$ (see figure 7b), instead of remaining nearly unchanged. This may be responsible for the optimal condition moving towards the lower bound of $O(1)$.

The energy output $C_p$ defined in (3.14) and presented in figure 7(a) is a quantity averaged over the length of the flag. Since the curvature of the fluttering flag is generally non-uniform, the generated energy varies along its length. To compare the contributions from different sections of the flag, we present in figure 9 the distributions of local curvature $\langle |\varLambda |\rangle$ and corresponding local energy output $\langle \upsilon ^2/\beta \rangle$ along the length of the flag, both being time-averaged quantities, for the case with $\beta =0.3$. Not surprisingly, they both monotonically increase with $s$ and reach their respective maximum at the trailing edge ($s=1$), consistent with what has been reported in Singh, Michelin & De Langre (2012) and Piñeirua, Doaré & Michelin (2015). The downstream half of the flag ($0.5< s<1$) contributes approximately $90\,\%$ of the overall energy, suggesting that in practice it is more cost effective to deploy piezoelectric patches only in this part.

Figure 9. Distributions of the energy output $\langle \upsilon ^2/\beta \rangle$ and the curvature $\langle |\varLambda |\rangle$ along the length of the flag for the case with $\beta =0.3$, where $\langle {\cdot }\rangle$ is the time-averaging operator.

4.2.2. The RLC circuit

When an inductor is introduced into the circuit, the circuit becomes an RLC circuit with an undamped natural frequency $\omega _0$. We have shown in § 4.1 that varying $\omega _0$ may result in mode transitions for the inverted piezohydroelastic flag, which is associated with the ‘lock-in’ phenomenon. In this subsection, we further evaluate its influence on the energy harvesting performance.

Figure 10 presents the variations of the flag energy output $C_p$ against the frequency ratio $\omega _0/\omega$ at different $\beta$ values. It is seen that, regardless the flag stiffness, the implementation of small $\omega _0$ (corresponding to very large inductance), e.g. $\omega _0/\omega =0.1$, makes the circuit close to a pure RC one, so that the energy output remains almost the same as that at $\omega _0=0$. On the other hand, very large $\omega _0$ tends to short circuit the loop, resulting in almost zero energy output. Similarly, when the resistance $\beta$ is too small or too large (e.g. $\beta =0.01$ or 100), the electrical loop is close to short circuited or a pure inductor-capacitor circuit. In either case the energy output is trivial. As such, the enhancement of energy output from the RLC circuit can only be possible when both $\omega _0$ and $\beta$ are moderate. Indeed, significant enhancement is observed for both the moderate flexible flag (i.e. $K^*=0.25$ in figure 10a) and the less flexible flag (i.e. $K^*=0.4$ in figure 10b), in the ‘lock-in’ regime where $\omega _0/\omega \sim O(1)$. For the moderate flexible flag, however, three exception cases are observed, i.e. $\beta =0.01$, 0.1 and 0.4, where the electrical circuit is overdamped and hence cannot induce larger structure motions through the ‘lock-in’. For the stiff flag, the performance improvement induced by the ‘lock-in’ seems more pronounced. This is attributed to the extremely low energy output produced in the small $\omega _0$ cases which are mainly in the unfavourable small-deflection mode (evidenced in figures 5 and 6).

Figure 10. Variation of $C_p$ against $\omega _0/\omega$ for (a) the moderate flexible flag ($K^*=0.25$) and (b) the less flexible flag ($K^*=0.4$). All peaks are marked with filled symbols. (c) Optimal frequency ratio $(\omega _0/\omega )_{opt}$ vs damping ratio $\zeta$. The scatters are $(\omega _0/\omega )_{opt}$ extracted from (a,b).

Upon closer inspection of figure 10(a,b), it is found that all maximum $C_p$ values (marked with filled symbols) occurring in the RLC-circuit regime (where $\omega _0 \neq 0$) are located slightly to the right of the $\omega _0/\omega =1$ line due to the damped nature of the circuit. The corresponding optimal frequency ratio $(\omega _0/\omega )_{opt}$ that leads to the maximum $C_p$ is expected to be highly associated with the circuit's damping ratio $\zeta =1/2\beta \omega _0$. As such, we map all these maximum $C_p$ values in the $(\omega _0/\omega )_{opt}-\zeta$ space as shown in figure 10(c). Two sequences of data labelled ‘a’ and ‘b’ are respectively extracted from panels (a) and (b). As expected, all these RLC-circuit boosted maximum $C_p$ values appear in the underdamped zone (i.e. $\zeta <1$). Moreover, all the corresponding optimal frequency ratios are greater than 1, and they generally grow with $\zeta$. In contrast, all the maximum $C_p$ values occurring in the RC-circuit regime (where $\omega _0 = 0$) are located at the same place in the overdamped zone (i.e. $\zeta >1$), where $\zeta =+\infty$ and $(\omega _0/\omega )_{opt}=0$.

To establish a better understanding of the optimal performance of the system, the DMD analysis is also conducted for the flag equipped with oscillatory circuits. We present in figure 11 the DMD analysis results of a representative sequence of cases with various $\omega _0$ and $\beta =1.0$ for the moderate flexible flag, which specifies the spectrum and energy output of the leading few DMD modes (calculated with (B11)), and compares the total energy output $C_p$ with the contributions from the first two DMD modes. As revealed in (a–e), when $\omega _0$ is small the energy output is dominated by the leading mode. However, as $\omega _0$ exceeds $1.0$ the energy output of the second modes (appearing at the third harmonic frequency) significantly increases, which can be even higher than the first mode (e.g. at $\omega _0=5.0$ in e). As such, in this $\omega _0$ regime the electrical dynamics is determined by the coupled effect of the first few modes, instead of solely by the leading mode. It is necessary to point out that, when higher frequency modes play a considerable role in the electrical dynamics, the linear combination of $C_{p,1}$ and $C_{p,2}$ may not reproduce $C_p$ exactly. This is because when there are multiple modes of the same order of magnitude contributing to the energy output, the interaction among these modes could be strong and hence cannot be neglected. Nevertheless, in the ‘lock-in’ regime where $\omega _0$ is near or less than 1.0, the leading DMD mode still plays a dominant role in contributing to $C_p$.

Figure 11. (a–e) The DMD spectra of the flag energy harvesting process and the corresponding energy output for different $\omega _0$ when $\beta =1.0$ and $K^*=0.25$. (f) Total energy output $C_p$ and contributions from the first or two DMD modes vs $\omega _0$. Panels show (a) $\omega _0=0$, (b) $\omega _0=0.5$, (c) $\omega _0=1.0$, (d) $\omega _0=2.0$, (e) $\omega _0=5.0$.

With the above observation, it is inferred that the $C_{p,1}$ formulation is still useful in determining the condition for optimal energy output. According to Appendix B, we have

(4.2)\begin{equation} C_{p,1}\sim\dfrac{\alpha^2K^*}{\beta\left(\dfrac{\omega_0^2}{{\rm Im}\{\omega_1\}^2}-1\right)^2+\dfrac{1}{\beta {\rm Im}\{\omega_1\}^2}}|b_1|^2. \end{equation}

As noted before, the amplitude $|b_1|$ and frequency $\textrm {Im}\{\omega _1\}$ of the leading DMD mode vary slightly in comparison with $\beta$ in the ‘lock-in’ regime. The denominator of the above formula is minimized approximately when

(4.3)\begin{equation} \beta\left(\dfrac{\omega_0^2}{{\rm Im}\{\omega_1\}^2}-1\right)^2 \approx \dfrac{1}{\beta {\rm Im}\{\omega_1\}^2}, \end{equation}

which can be rewritten as $\textrm {Im}\{\omega _1\}^4-(2\omega _0^2+1/\beta ^2)\textrm {Im}\{\omega _1\}^2+\omega _0^4\approx 0$. Solving it gives

(4.4)\begin{equation} {\rm Im}\{\omega_1\}^2\approx\omega_0^2\left(1+2\zeta^2\pm2\zeta\sqrt{1+\zeta^2} \right), \end{equation}

where $\zeta =1/(2\beta \omega _0)$ is the damping ratio of the electrical circuit. We obtain a simple expression of the optimal frequency ratio that leads to the maximum energy output

(4.5)\begin{equation} \left(\dfrac{\omega_0}{\omega}\right)_{opt}\approx\dfrac{1}{\sqrt{1+2\zeta^2\pm2\zeta\sqrt{1+\zeta^2}}}=\varPi(\zeta). \end{equation}

If plotting (4.5) by adopting the negative sign in figure 10(c), we see all the optimal frequency ratios appearing in the RLC-circuit regime (see figure 10a,b) are located near the prediction (the black solid curve) in the underdamped zone, following a similar trend. Due to the damping effect exerted by the circuit on the flag dynamics, all these values are located above the prediction curve. It is interesting to see that (4.5), if adopting the positive sign, also captures all the optimal frequency ratios appearing in the RC-circuit regime, even though they collapse at the same point in the overdamped zone. Overall speaking, (4.5) provides a fairly good prediction of the optimal frequency for the electrohydroelastic flag where the electrical damping ratio is known.

Although very insightful, (4.5) includes the system response frequency $\omega$, which is unknown beforehand. To better facilitate the selection of electrical circuit for energy harvesting, especially when the electrical load ($\beta$) is known and the inductance ($\ell$ or $\omega _0$) is to be determined, $C_p$ contours for the moderate and less flexible flags are plotted in the $\omega _0-\beta$ space, as shown in figure 12(a,b), respectively, in which all the optimal cases (with maximum $C_p$ for given $\beta$) are marked with red circles. To underline the potential of performance enhancement by using RLC circuits, figure 12(c,d) further compare the maximum $C_p$ values at various $\beta$ when the RC or RLC circuit is deployed, which are extracted from the contours. Obviously, at lower $\beta$ (the range varies with the stiffness of the flag), the maximum $C_p$ is achieved with RC circuits, indicating that the employment of an inductor fails to improve the system performance. At higher $\beta$, on the contrary, the best performance is achieved by RLC circuits (i.e. with non-zero $\omega _0$) through the ‘lock-in’. The boundary that separates these two optimal regimes, as shown in the figure, is the $\zeta =1.0$ (the critical damping) line. The implication is that, to trigger the ‘lock-in’ state, one should carefully choose the inductance $\omega _0$ to ensure the damping ratio falls into the underdamped zone, i.e. $\omega _0>1/2\beta$, which defines the lower bound of the optimal $\omega _0$.

Figure 12. (a,b) Contour map of $C_p$ in the $\beta -\omega _0$ space and (c,d) comparison of $C_p$ between the optimal RLC circuit (extracted from the contour maps along the red trajectories) and the RC circuit: (a,c) the moderate flexible flag and (b,d) the less flexible flag. Panels show (a) $K^*=0.25$, (b) $K^*=0.4$, (c) $K^*=0.25$, (d) $K^*=0.4$.

It is also observed from figure 12 that, in the underdamped zone (i.e. $\zeta <1$), the maximum $C_p$ generally decreases with $\beta$ despite some local peak or plateau. This suggests a possible upper bound for the optimal $\omega _0$. From (4.5), we can obtain in the underdamped zone the maximum value of optimal $\omega _0/\omega$, i.e. $\max (\omega _0/\omega )_{opt}\approx 2.4$, which corresponds to the intersection point value of the (4.5) curve and the $\zeta =1$ line, as denoted in figure 10(c). Note that, when undergoing the flutter modes, the flapping frequency of the electrohydroelastic flag $\omega$ is upper bounded by the frequency of the pure mechanical flag $\omega _s$, i.e. $\omega <\omega _{s}$, which is evidenced in figure 6(a,c). This leads to a reasonable upper bound for optimal $\omega _0$, i.e. $\omega _0<2.4\omega _s$. Therefore, the optimal $\omega _0$ shall be chosen in a double-bounded range as

(4.6)\begin{equation} 0.5\beta^{-1}<\omega_0<2.4\omega_{s}. \end{equation}

This bounded range is depicted by the triangle zone formed with the two black dashed lines shown in figure 12(a,b). Note that some optimal $\omega _0$ value slightly falls outside this triangle zone, especially near the $\zeta =1$ line (see figure 12b). This can be attributed to the general underestimation made by (4.5) through ignoring the damping the circuit brings to the electrohydroelastic flag. Nevertheless, the bounded range described by (4.6) still delineates an acceptable zone for optimal $\omega _0$.

Note that (4.6) requires a condition

(4.7)\begin{equation} \beta > 0.21\omega_s^{-1}. \end{equation}

It indicates that, if one would like to utilize RLC circuits to achieve better performance over RC circuits, the electrical load $\beta$ must be greater than a critical value, i.e. $0.21\omega _s^{-1}$. Otherwise, the use of RLC circuits cannot bring any benefit.

The optimal frequency ratio (4.5), the proposed searching zone for the optimal $\omega _0$ (4.6) and the critical electrical load (4.7) are very useful in guiding the optimization process of the electrical circuits. They can assist in choosing an appropriate inductor quickly and efficiently when RLC circuits are needed for higher energy production, which usually happens at relatively large electrical load.

It should be pointed out that, although the free-stream velocity $U^*$ is not explicitly adopted as a design variable in the above optimization process, its influence has already been reflected in (4.5)–(4.7). Specifically, in (4.5) it affects $\omega$, the actual response frequency of the system, through complex fluid–structure–electric coupling, whereas in (4.6) and (4.7) it affects the ranges of $\omega _0$ and $\beta$ through $\omega _s$, the flutter frequency of the pure mechanical flag.

5. Conclusions and discussion

The present study focuses on the multi-physics dynamics of an inverted electrohydroelastic flag. The mechanical energy of ambient flows is harvested through FIV, which is then converted into electricity via an electric circuit. By using a fully coupled fluid–structure–electric model, we simulate the entire energy conversion process, with the focus placed on the dynamics of the flag and its energy output, particularly the impacts of structural and electrical parameters.

Via parametric studies, a variety of structural response modes are identified. Among these response modes, the two flutter modes, i.e. the symmetric- and asymmetric-flutter modes, are the most suitable ones for sustainable energy harvesting. We found that, with moderate flexibility, the flag is more likely to fall into the flutter modes, whereas too flexible flags can be easily bent and deflected to one side considerably (large-deflection mode) and too rigid flags tend to remain almost undeformed (small-deflection mode). We also found that, with only an RC circuit deployed, the increase of electrical load (i.e. $\beta$) can stabilize the flag.

For benchmark purposes, the energy harvesting performance of the flag equipped with an RC circuit is first examined. Our simulation results and DMD-analysis-based derivations confirm the optimal condition, i.e. $\beta \omega \sim O(1)$, for the maximum energy output. However, unlike in the weak electromechanical coupling scenario, where this optimal condition works very well, the strong coupling considered in the present study connotes more complicated nonlinear fluid–structure–electric interaction and pushes this condition towards the lower bound of $O(1)$.

We then proceed to study the impact of inductors on the performance of the system. By adding an inductor, the output loop becomes an oscillatory RLC circuit with an undamped natural frequency $\omega _0$. Our results show that, when the circuit's natural frequency $\omega _0$ satisfies $0.5\beta ^{-1}<\omega _0<2.4\omega _{s}$ (i.e. (4.6)), the employment of the oscillatory circuit tends to lock the structural vibration frequency into the circuit's frequency, which can significantly improve the energy harvesting performance. This performance improvement due to the ‘lock-in’ can only occur in the underdamped zone. Furthermore, the optimal $\omega _0$ can be estimated by (4.5). These findings and guidelines are expected to play an important role in the RLC circuit design and optimization.

In the present study, an RLC circuit is used to harvest the energy from ambient flows. However, we note that the simplified electrical rendition is still far from the real system to be developed. For example, for voltage compatibility, electrical interfaces need to be applied between the terminal electric load and the piezoelectric element (Sarker et al. 2019). The most frequently used interface is the combination of a diode rectifier bridge and a filter capacitor, which are respectively used for rectifying and smoothing the alternating voltage (Lefeuvre et al. 2006). Moreover, more complicated interfaces based on the so-called synchronized switching harvesting with inductor technique based on a specific nonlinear processing of the voltage have been proposed (Daniel et al. 2016; Li, Roy & Calhoun 2019). Such interfaces can significantly affect the dynamics of the flag and influence the energy output. Therefore, a complete multi-physics model that involves an electric interface may be a better option, whose results can be directly used for guiding the prototype development. This will be our future direction in this area.