2.1 Configuration

As shown in figure 1, we investigate flow-induced vibration of an inverted piezoelectric flag with length $L$ , placed in a uniform steady flow with velocity $U$ and density ${\it\rho}$ , with the initial inclination angle to the flow of ${\it\psi}$ . The flag is assumed to be inextensible, free at its leading edge and clamped at the trailing edge. It is mathematically represented as a nonlinear 2D Euler–Bernoulli beam and the position of its midline is denoted with $\boldsymbol{X}(s,t)$ , where $t$ is the time and $s$ is the Lagrangian coordinate along the beam, starting from $s=0$ at the leading edge and ending at $s=L$ at the trailing edge.

Figure 1.

Schematic configuration of the inverted flag in uniform flow.

2.2 Governing equations for flow–structure interaction

The flow is governed by the incompressible Navier–Stokes equations, which in non-dimensional form are expressed as

where the density of the fluid $({\it\rho})$ , the length of the flag $(L)$ and the free-stream velocity $(U)$ are respectively chosen to be the characteristic density, length and velocity in this study and $Re={\it\rho}UL/{\it\mu}$ is the Reynolds number. In the above equation, $\boldsymbol{x}$ is the location vector in the fluid domain, $\boldsymbol{u}$ is the fluid velocity, $p$ is the dynamic pressure, ${\it\rho}$ is the density of the fluid, ${\it\mu}$ is the dynamic viscosity and $\boldsymbol{f}$ is the momentum forcing term associated with boundary conditions on the surface of the flag.

The flag is chosen to be inextensible, with the inextensibility constraint expressed as $\text{d}(|\partial X/\partial s|)/\text{d}t=0$ (Tornberg & Shelley Reference Tornberg and Shelley2004; Shoele & Zhu Reference Shoele and Zhu2012), resulting in the following system of equations for the position and tension force of the flag:

Here ${\bf\tau}=\partial \boldsymbol{X}/\partial s$ is the unit tangent vector along the axial coordinate $s$ , $\boldsymbol{n}$ is the unit normal vector along the axial coordinate, and $m_{s}$ is the excess mass per unit length of the flag defined as $m_{s}=m_{f}-{\it\rho}A_{s}$ , where $m_{f}$ is the mass per unit length of the flag, ${\it\rho}$ is the density of the fluid and $A_{s}$ denotes the cross-sectional area of the flag. In the limiting case of $A_{s}\rightarrow 0$ (a very thin flag), which is relevant to the current configuration, $m_{s}\sim m_{f}$ . Furthermore, ${\it\sigma}$ is the tension, ${\it\gamma}$ is the transverse stress and $\boldsymbol{F}$ is the Lagrangian forcing exerted on the flag by the surrounding fluid. The free- and clamped-end boundary conditions are imposed at the leading and trailing edges of the flag, respectively, as

where ${\it\psi}$ is the inclination of the unperturbed flag to the flow direction.

2.3 Governing equations for piezo–structure coupling

Following Michelin & Doare (Reference Michelin and Doare2013), we assume that both surfaces of the flag are covered with infinitesimal piezoelectric patches with segmentation lengths much smaller than $L$ . Moreover, it is assumed that the output electric circuit is coupled to the deformation of the flag through mutual interaction between the electric voltage of the patches, local stretch or compression of the patches as a result of the beam deflection and charge transfer through the thickness of the flag. Considering the limiting case of continuous coverage by the piezoelectric patches, we can further express the coupling between the electric voltage between the positive electrodes of each infinitesimal patch $V(s,t)$ , the charge transfer per unit length $Q$ and the flag deformation as (Michelin & Doare Reference Michelin and Doare2013)

where $c$ is the linear capacitance of the piezoelectric element, and ${\it\chi}$ is the coupling coefficient related to the intrinsic material properties of piezoelectric patches. In this study, we assume a purely resistive circuit between $Q$ and $V$ , for which the electrical equation of state can be written as

where ${\it\varsigma}$ is the linear conductivity coefficient between the piezoelectric patches on the top and bottom surfaces of the flag, which is equal to the inverse of the effective resistance of a purely resistive circuit connecting the top and bottom patches. As per (Michelin & Doare Reference Michelin and Doare2013) the total moment $M$ along the flag can be expressed as

where $k_{b}$ is the bending rigidity. The Euler–Bernoulli assumption (Audoly & Pomeau Reference Audoly and Pomeau2010) for a slender piezoelectric flag leads to

Denoting the characteristic voltage and charge density with $U\sqrt{{\it\rho}L/c}$ and $U\sqrt{{\it\rho}Lc}$ , respectively, we can write (2.3)–(2.10) in their non-dimensional forms. In the following, we use the same notations for the dimensionless quantities as their dimensional counterparts and further combine (2.10) and (2.3) to reach the final non-dimensional equations for the piezoelectric flag as

where ${\it\sigma}$ is also calculated from the non-dimensional counterpart of (2.4). Five non-dimensional parameters characterize the coupled fluid–structure–electrical system. They are the Reynolds number $Re$ , the inertial ratio between the fluid and structure $M^{\ast }$ , the non-dimensional free-stream velocity $U^{\ast }$ , the coupling coefficient ${\it\alpha}$ , and the tuning coefficient of the electrical system ${\it\beta}$ , which are defined as

2.4 Immersed boundary method

Following the conventional approach to the immersed boundary method (Peskin Reference Peskin1972), the Lagrangian velocity $\boldsymbol{U}(s,t)$ and force density $\boldsymbol{F}(s,t)$ at the position $s$ along the flag are related to $\boldsymbol{u}(\boldsymbol{x},t)$ and $\boldsymbol{f}(\boldsymbol{x},t)$ through

respectively, where ${\it\Omega}_{F}$ is the fluid domain, ${\it\Gamma}$ the surface of the flag and ${\it\delta}$ the Dirac delta function.

In the present immersed boundary method, the fluid and piezoelectric flag equations are solved independently, and their interaction force is formulated in terms of the penalty momentum forcing (Huang & Sung Reference Huang and Sung2009). Two groups of points are considered at the place of the flag: one group is the actual flag Lagrangian points, which are assumed not to interact with the fluid directly; the other group is massless virtual points that are assumed to move with the local fluid velocity at their position. The two groups of points are connected together with stiff springs and dampers, to model the interaction between the structure and the fluid, and therefore the momentum forcing $\boldsymbol{F}(s,t)$ is expressed as

where ${\it\kappa}_{p}$ is the large penalty constant, ${\rm\Delta}t$ is the computational time step, $\tilde{\boldsymbol{X}}$ and $\tilde{\boldsymbol{U}}$ are the position and velocity of massless virtual points, respectively, and $\boldsymbol{X}$ and $\boldsymbol{U}$ are the position and velocity of the Lagrangian points on the flag. The position $\tilde{\boldsymbol{X}}$ and velocity $\tilde{\boldsymbol{U}}$ are computed from the current and previous fluid velocities as

This formulation is identical to what has been proposed by Goldstein, Handler & Sirovich (Reference Goldstein, Handler and Sirovich1993) as the feedback forcing scheme, and has successfully been used in immersed boundary methods for a variety of applications (Huang, Shin & Sung Reference Huang, Shin and Sung2007; Huang & Sung Reference Huang and Sung2009; Shoele & Zhu Reference Shoele and Zhu2012, Reference Shoele and Zhu2013; Shoele & Mittal Reference Shoele and Mittal2014). In particular, in the feedback forcing scheme, we have

and it can be seen that (2.16) is indeed a special case of the feedback forcing method with $a_{1}={\it\kappa}_{p}$ and $a_{2}={\it\kappa}_{p}{\rm\Delta}t$ , with a benefit that, in the current method, there is only one penalty parameter, which takes care of both the velocity and position effects.

The formulations in (2.11) and (2.12) are solved using the finite-difference approach. The fluid equations are solved on a staggered grid, where the velocity and momentum forcing are defined on the staggered grid and the pressure is defined at the centres of cells. Further details of this formulation can be found in Shoele & Zhu (Reference Shoele and Zhu2012).

2.5 Numerical implementation

1. (i) Fluid equations:

   (2.20) $$\begin{eqnarray}\displaystyle & & \displaystyle \frac{\boldsymbol{u}^{n+1,k+1}-\boldsymbol{u}^{n}}{{\rm\Delta}t}+\frac{1}{2}(N(\boldsymbol{u}^{n+1,k+1})+N(\boldsymbol{u}^{n}))\nonumber\\ \displaystyle & & \displaystyle \quad =-\boldsymbol{{\rm\nabla}}p^{n+1/2,k+1}+\frac{1}{2Re}{\rm\nabla}^{2}(\boldsymbol{u}^{n+1,k+1}+\boldsymbol{u}^{n})+\boldsymbol{f}^{n+1,k},\end{eqnarray}$$
   (2.21) $$\begin{eqnarray}\displaystyle \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}^{n+1,k+1}=0. & & \displaystyle\end{eqnarray}$$

2. (ii) Structural equations:

   (2.22) $$\begin{eqnarray}\displaystyle \frac{\boldsymbol{X}^{n+1,k+1}-\boldsymbol{X}^{\ast }}{{\rm\Delta}t^{2}} & = & \displaystyle D_{s}({\it\sigma}^{n+1/2,k+1}D_{s}\boldsymbol{X}^{n+1,k+1})-D_{ss}\left(\frac{M^{\ast }}{U^{\ast ^{2}}}D_{ss}\boldsymbol{X}^{n+1,k+1}\right)\nonumber\\ \displaystyle & & \displaystyle +\,D_{s}\left(\frac{{\it\alpha}\sqrt{M^{\ast }}}{U^{\ast }}D_{s}v^{n+1,k}\boldsymbol{n}^{n+1,k}\right)+M^{\ast }\boldsymbol{F}^{n+1,k},\end{eqnarray}$$
   (2.23) $$\begin{eqnarray}\displaystyle D_{ss}({\it\sigma}^{n+1/2,k+1}{\bf\tau}^{n+1,k})\boldsymbol{\cdot }{\bf\tau}^{n+1,k} & = & \displaystyle \frac{1-2({\bf\tau}\boldsymbol{\cdot }{\bf\tau})^{n}+({\bf\tau}\boldsymbol{\cdot }{\bf\tau})^{n-1}}{2{\rm\Delta}t^{2}}-(D_{s}\boldsymbol{U}\boldsymbol{\cdot }D_{s}\boldsymbol{U})^{n+1,k}\nonumber\\ \displaystyle & & \displaystyle -\,D_{s}\left(D_{s}\left(\frac{M^{\ast }}{U^{\ast ^{2}}}{\it\gamma}\boldsymbol{n}\right)^{n+1,k}-M^{\ast }\boldsymbol{F}^{n+1,k}\right)\boldsymbol{\cdot }{\bf\tau}^{n+1,k},\nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$
   (2.24) $$\begin{eqnarray}\displaystyle \boldsymbol{U}^{n,k+1}=\frac{\boldsymbol{X}^{n+1,k+1}-\boldsymbol{X}^{n}}{{\rm\Delta}t}. & & \displaystyle\end{eqnarray}$$

3. (iii) Electrical equation:

   (2.25) $$\begin{eqnarray}\displaystyle {\it\beta}\frac{V^{n+1,k+1}-V^{n}}{{\rm\Delta}t}=-V^{n+1,k+1}-\frac{{\it\alpha}{\it\beta}\sqrt{M^{\ast }}}{U^{\ast }}\frac{D_{ss}\boldsymbol{X}^{n+1,k+1}\boldsymbol{\cdot }\boldsymbol{n}^{n+1,k}-D_{ss}\boldsymbol{X}^{n}\boldsymbol{\cdot }\boldsymbol{n}^{n}}{{\rm\Delta}t}.\qquad & & \displaystyle\end{eqnarray}$$

4. (iv) Fluid–structure interaction:

   (2.26) $$\begin{eqnarray}\displaystyle & \displaystyle \tilde{\boldsymbol{U}}^{n+1,k+1}=\mathop{\sum }_{x}\boldsymbol{u}^{n+1,k+1}(\boldsymbol{x}){\it\delta}_{h}(\boldsymbol{x}-\tilde{\boldsymbol{X}}^{n+1,k}){\rm\Delta}x{\rm\Delta}y, & \displaystyle\end{eqnarray}$$

   (2.27) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\tilde{\boldsymbol{X}}^{n+1,k+1}-\tilde{\boldsymbol{X}}^{n}}{{\rm\Delta}t}=\tilde{\boldsymbol{U}}^{n+1,k+1}, & \displaystyle\end{eqnarray}$$

   (2.28) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{F}^{n+1,k+1}=-{\it\kappa}_{p}\frac{{\it\rho}U^{2}}{L}[(\tilde{\boldsymbol{X}}^{n+1,k+1}-\boldsymbol{X}^{n+1,k+1})+{\rm\Delta}t(\tilde{\boldsymbol{U}}^{n+1,k+1}-\boldsymbol{U}^{n+1,k+1})], & \displaystyle\end{eqnarray}$$

   (2.29) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{f}^{n+1,k+1}=\mathop{\sum }_{s}\boldsymbol{F}^{n+1,k+1}(s){\it\delta}_{h}(\boldsymbol{x}-\boldsymbol{X}^{n+1,k+1}){\rm\Delta}s. & \displaystyle\end{eqnarray}$$
   In the above equations, $N$ is the discrete convective operator, $D_{s}$ is the difference operator with respect to $s$ , ${\rm\Delta}x$ and ${\rm\Delta}y$ are the sizes of the uniform fluidic grid around the flag, ${\rm\Delta}s$ is the size of the Lagrangian grid along the flag, and ${\it\delta}_{h}$ is chosen to be a four-point smoothed delta function (Peskin Reference Peskin2002),
   (2.30) $$\begin{eqnarray}\displaystyle {\it\delta}_{h}(\boldsymbol{x})=\frac{1}{{\rm\Delta}x{\rm\Delta}y}\,{\it\varphi}\left(\frac{x}{{\rm\Delta}x}\right){\it\varphi}\left(\frac{y}{{\rm\Delta}y}\right), & & \displaystyle\end{eqnarray}$$
   with
   (2.31) $$\begin{eqnarray}\displaystyle {\it\varphi}(r)=\left\{\begin{array}{@{}ll@{}}\displaystyle (1/8)\left(3-2|r|+\sqrt{1+4|r|-4r^{2}}\right), & 0\leqslant |r|<1,\\ \displaystyle (1/8)\left(5-2|r|+\sqrt{-7+12|r|-4r^{2}}\right), & 1\leqslant |r|<2,\\ 0, & 2\leqslant |r|.\end{array}\right. & & \displaystyle\end{eqnarray}$$

We use a second-order approximation of the nonlinear term $N$ in (2.20), obtained from the approximate second-order temporal factorization of (2.1). The resultant systems of tri-diagonal matrices for the intermediate velocities are then solved using the efficient Thomas algorithm in alternative directions. The velocities and pressure are decoupled by using the block LU decomposition (Kim, Baek & Sung Reference Kim, Baek and Sung2002), and the pressure Poisson equation is solved with the BiCGSTAB iterative method along with the ILU $(k)$ preconditioner (Saad Reference Saad1990). For the initialization of the fixed-point iteration at $k=0$ for a given time step, we use the results from the previous time step. The iterative procedure continues until the values are converged down to $O(10^{-8})$ .

Figure 2.

The computational domain and assumed boundary conditions.

Figure 2 shows the rectangular computational domain as well as the specified boundary conditions employed in this study. The initial position of the flag is a straight line extending from its leading edge to its trailing edge at $(L,0)$ . The flag may be inclined at an angle ${\it\psi}$ to the incoming flow but the value of ${\it\psi}$ is equal to $0^{\circ }$ unless noted otherwise. The majority of the simulations are conducted at a Reynolds number of $Re=200$ . This value is large enough to generate complex vortex-induced vibrations in the flag, but low enough to allow for a thorough sweep of the parameter space. However, one part of the study also examines the effect of Reynolds number on the configuration. A uniform velocity is assumed at all the boundaries except the outflow boundary, where a convective outflow boundary condition is used. A zero normal gradient boundary condition is used for pressure at all the boundaries.

The size of the computational domain is chosen to be $-8L\leqslant x\leqslant 27L$ and $-10L\leqslant y\leqslant 10L$ after a domain dependence study to ensure that further increase in domain size has no significant effect on the gross quantities (lift, drag, vibration amplitude and frequency) as well as the vortex topology in the vicinity of the flag. A uniform grid with 78 points is used to discretize the flag and the computational domain is discretized with $555\times 420$ points in the $x$ and $y$ directions, respectively. The grid density is increased in the region around the flag extending from $-L$ to $2L$ in the $x$ direction and $-1.5L$ to $1.5L$ in the $y$ direction (as shown by the a shaded area in figure 2) with grid sizes of ${\rm\Delta}x={\rm\Delta}y=0.0167L$ . The time step $t$ is chosen to be $0.00025L/U$ , which corresponds to a maximum Courant–Friedrichs–Lewy (CFL) number of approximately 0.04. In the numerical simulations we choose ${\it\kappa}_{p}=-2\times 10^{5}$ . Sensitivity tests conducted for the current system indicate that the chosen numerical parameters are adequate for the accuracy and stability of the simulations, and the maximum uncertainty introduced by these parameters is estimated to be less than 2 %.

2.6 Validation and verification

The accuracy and fidelity of the above formulation for flow and fluid–structure interaction has been assessed for a wide variety of canonical problems in the past (Shoele & Zhu Reference Shoele and Zhu2012, Reference Shoele and Zhu2013, Reference Shoele and Zhu2015; Shoele & Mittal Reference Shoele and Mittal2014). In particular, previous simulations of flow past a flexible flag have shown that the immersed boundary solver employed here accurately predicts vortex-dominated flows at the Reynolds numbers relevant to the current study. One verification study that is particularly relevant to the current work is the one in Shoele & Mittal (Reference Shoele and Mittal2014), where the flow-induced flutter of a flexible plate mounted behind a circular cylinder in a channel was simulated and the results compared to a study that employed a finite-element discretization (Turek & Hron Reference Turek and Hron2006) and another study that employed a sharp-interface immersed boundary method (Bhardwaj & Mittal Reference Bhardwaj and Mittal2012). The comparisons showed that the results generated by the current immersed boundary flow–structure interaction (FSI) solver were in very good agreement with these past studies, thereby verifying the fidelity of the current solver for an FSI problem similar to the one studied here. The details of these studies are provided in the references quoted above and will not be repeated here.

The new modelling component that is incorporated in the current study is the structural–electrical coupling of the piezoelectric beam. In order to test this module of the computational model, we employ the canonical cases of a thin piezoelectric plate with one end cantilevered, undergoing forced oscillations in the absence of an ambient fluid. The plate is actuated with a periodic heaving motion at its fixed end of the form $y_{0}=A_{0}\sin ({\it\omega}_{0}t)$ and the vertical deflection of its tip examined. Under the assumption of small deflection, where nonlinear effects can be neglected, and the base of the beam actuated with a periodic function, the vertical deflection of the beam can be written as $y(x,t)=\text{Re}[Y(x)\text{e}^{\text{i}{\it\omega}t}]$ , and we can combine (2.11) and (2.12) into

where

along with the boundary conditions at the free and cantilevered ends, which are $Y_{xxx}(1)=Y_{xx}(1)=Y_{x}(0)=0$ and $Y(0)=1$ . The ordinary differential equation (2.32) is solved analytically and the deflection of the beam at its tip is found to be

The above then gives us an analytical solution of the coupled mechanical–electrical problem in the case of small deflection for a simple canonical case and we can obtain solutions for any choice of the electrical and structural parameters. The same problem can also be solved numerically by using our computational model.

In figure 3, we compare our prediction of the fundamental frequency of the beam (with length $L$ ) as well as the amplification factor of the tip oscillation with the analytical results for a configuration with $A_{0}/L=0.0025$ , $M^{\ast }=1$ and $U^{\ast }=1$ , but with two different electrical properties: one with ${\it\alpha}=0.3$ and ${\it\beta}=1$ , and the other with ${\it\alpha}=1.0$ and ${\it\beta}=0.3$ . The comparison shows that the computational method captures the overall trends quite well and does a reasonable job of predicting the resonance frequencies of the coupled piezo–structural system. A trend that emerges from the comparison is that the peaks from the simulations are shifted down and slightly to the right of those predicted by the analytical model. Further examination of this shift in the computational predictions is due to nonlinear effects. Note that even though the input deflection is small, the deflection in the beam at resonance can be amplified by factors ranging from approximately 10 to 200, and these are large enough to interject nonlinear effects into the dynamics of this beam; these nonlinear effects are absent in the analytical solution but are present in the computational model, leading to the observed differences. Stanton et al. (Reference Stanton, Erturk, Mann, Dowell and Inman2012) showed that, for large deflections, the nonlinearity of the mechanical–electrical coupling term reduces the response amplitude at resonance and the inextensibility condition of the beam introduces a cubic nonlinearity in the bending stiffness, which shifts the response curves to the left. Thus, the above comparison provides some confidence in the fidelity of the piezo–structural coupling and also provides a preview of the non-negligible effect of this coupling on the dynamics of the flexible piezoelectric beam.

Figure 3.

Spectrum of the amplification factor at the tip of the piezoelectric beam $(A_{tip}/A_{base})$ as a function of actuation frequency at the base of the beam. The dotted lines are the theoretical prediction for a small-deflection model and the solid line along with the symbols are the numerical results.

2.7 Key metrics for flag dynamics and performance

A variety of metrics are employed in order to characterize the dynamics of the flag as well as its energy harvesting performance. The mean deflection of the flag is defined by $\bar{{\it\theta}}=\tan ^{-1}(\bar{y}_{LE}/\bar{x}_{LE})$ , where $x_{LE}$ and $y_{LE}$ are the $x$ and $y$ locations of the leading edge of the inverted flag, respectively, and the overbar indicates a time average. The non-dimensional oscillation amplitude $A_{max}$ of the flag leading edge is defined as

where ‘primed’ values denote root-mean-square (r.m.s.). The mean value and r.m.s. are calculated over a period of approximately four vibration cycles once the system reaches a stationary state.

The energy harvested by the piezoelectric flag is equal to the instantaneous power dissipated in the piezo patches. Following Michelin & Doare (Reference Michelin and Doare2013), the harvested energy is calculated in terms of a power coefficient by normalizing the harvested power with ${\it\rho}U^{3}(2L)/2$ , which is the kinetic energy flux of the fluid in the region accessible to the flag:

This power coefficient can also be viewed as the aerodynamic efficiency of energy harvesting by the flag.