Electro-Aero-Mechanical Model of Piezoelectric Direct-Driven Flapping-Wing Actuator

Abstract

We present an analytical model of a flapping-wing actuator, including its electrical, aerodynamic, and mechanical systems, for estimating the lift force from the input electrical power. The actuator is modeled as a two-degree-of-freedom kinematic system with semi-empirical quasi-steady aerodynamic forces and the electromechanical effect of piezoelectricity. We fabricated actuators of two different scales with wing lengths of 17.0 and 32.4 mm and measured their performances in terms of the stroke/pitching angle, average lift force, and average consumed power. The experimental results were in good agreement with the analytical calculation for both types of actuators; the errors in the evaluated characteristics were less than 30%. The results indicated that the analytical model well simulates the actual prototypes.

1. Introduction

Micro aerial vehicles (MAVs) have recently attracted significant attention owing to the increasing demand for automated transport, precision agriculture, and unmanned exploration/rescue missions. Among various types of MAVs, bioinspired flapping-wing robots have attracted particular interest owing to their potential to mimic the high agility and robustness of natural insects or birds. In particular, insect or hummingbird life-scale MAVs are considered useful because their extremely low weight and small size can enable safe flight and widespread applications.

There are many candidate actuation principles for flapping-wing MAVs, including piezoelectric, electromagnetic, shape memory alloy, and electrostatic. Piezoelectric actuation appears most suitable for small-scale vehicles owing to its high power-to-weight ratio, and some studies have already used piezoelectric actuators. For example, Harvard University [1,2,3] developed “RoboBee” and demonstrated stable tethered flight of a piezoelectric-actuated vehicle with a 30-mm wing span. This robot is equipped with piezoelectric Pb(Zr,Ti)O₃ (PZT) bimorph actuators and lever mechanisms to amplify and transmit the output displacement of the actuator to wings.

To realize fully untethered flight, an energy-efficient design is important. Because a flapping-wing actuator includes an electrical system, aerodynamic effects, and mechanical parts, a model integrating these three physics is required for detailed characterization. Several studies have reported mechanical-aerodynamic coupling models with nonlinear aerodynamic force [4,5] using semi-empirical quasi-steady aerodynamic theory [6,7,8,9]. The reported results showed that the models were in good agreement with experiments. In contrast to aerodynamics, the integration of electrical models has received less attention. For instance, a research group at UC Berkeley first presented an integrated model of a flapping-wing actuator that integrated the electric, aerodynamic, and mechanical parts to estimate the lift force and power consumption [10,11,12]. However, in their model, the aerodynamic force was approximated as a linear damping although it is basically nonlinear. The Harvard University research group considered the scaling of their robot on lift force, weight, and power consumption [13]. Because they focused on refining the design based on a specific prototype, their approach was heuristic, and they have not presented a deterministic model. Another group reported system modelling based on the measured response of the actuator [14]. This approach was very effective for designing a control system with an existing actuator; however, it was difficult to apply to the prediction of a newly designed actuator. A research group at the University of Southern California developed an MAV called “Robo Raven”, a larger, bird-scale vehicle, and reported the dynamic modelling of motor-driven flapping-wing MAVs and validation with flight data [15,16]. Another group at the California Institute of Technology developed “Robobat”, a bat-inspired MAV, and reported aerodynamics-based model optimization of the lift force [17]. However, these studies did not confirm the scalability of their models down to insect-scale vehicles. We believe that constructing an integrated model for an insect-scale actuation system can significantly contribute to the development of flapping-wing MAVs. Whitney et al. analyzed the conceptual design of insect-scale flapping-wing MAVs with a generalized (unspecified) actuator and aerodynamic force models [18]. They theoretically reported on how flight endurance was limited by power consumption without any experiments.

In light of the abovementioned previous studies, in the present study, we focused on the experimental validation of a flapping-wing actuator based on a piezoelectric actuation system. We present a model that integrates a flapping-wing actuator with a piezoelectric unimorph actuator and a several-tens-of-millimeter-long wings, as reported previously [19,20]. This actuator uses a direct-driven mechanism. In contrast to the mechanism used in RoboBee, this actuator does not have any amplifying transmission system; instead, a wing is directly connected to the unimorph actuator. This direct-driven mechanism was first proposed by a research group at the Army Research Laboratory [21,22], and it is advantageous in terms of structural simplicity. Two different-sized prototypes were fabricated and evaluated. We verified the theoretical model by comparing it to experimental results.

2. Electro-Aero-Mechanical Model

2.1. Direct-Driven Flapping-Wing Actuator

Figure 1 shows the actuation principle of a piezoelectric unimorph actuator and the structure and motion of the flapping-wing actuation system with the direct-driven mechanism. A unimorph actuator consists of a combination of a piezoelectric plate and a shim plate (elastic material), as shown in Figure 1a. Upon applying a voltage to the piezoelectric material, the actuator bends owing to the induced transverse piezoelectric stress. Figure 1b shows the structure of the flapping-wing actuation system. As described above, this system comprises only a unimorph piezoelectric actuator and a wing, and the wing is directly attached to the tip of the unimorph actuator. By actuating the unimorph, the wing is swung around the y-axis. This y-axis motion is simply called the “stroke motion”, and the stroke angle is expressed by $\varphi$. When the voltage signal frequency matches the resonance frequency of the system, large stroke amplitude is obtained. We define the x′y′z′-coordinate frame that rotates with the stroke motion. The wing is passively rotated around the x′-axis by air pressure associated with the stroke motion. This x′-axis rotation is called “pitching motion”, and its rotation angle is defined by $\psi$. We define the x′y″z″-coordinate as a wing fixed frame. The pitching motion is crucial to generate lift force; the pitched wing gets a certain angle of attack against the stroke motion and it generates lift force along the negative y-axis. Therefore, the wing should be designed to smoothly rotate around the x′-axis and to be rigid against the other directions. For this purpose, hinge structures have often been equipped near the top edge of the wing.

In the section below, we describe the analytical model of the actuation system. To derive this model, we employ several approximations for simplicity:

• The system has only two degrees of freedom (DOFs), $\varphi$ and $\psi$. Specifically, the wing can not only rotate but also translate on the x-z plane, and therefore, this assumption is generally incorrect. However, under a cyclic resonant condition, the wing motion can be approximated as a simple rotation around a certain axis. This approximation is verified based on the experimental results reported in a later section.
• Aerodynamics forces/moments are modeled as the sum of the translational force, force due to the added-mass effect, and rotational damping, and we use the blade element theory proposed in [5,6,7,8,9].
• We ignore the piezoelectric and other material losses and anchor-loss because their contribution is much smaller than that of aerial damping.

We derive the model based on Lagrange’s motion equation [23]. The total kinematic energy $T$ of the system is expressed as

$$T=T_A+T_B+T_w,$$

(1)

where $T_A$, $T_B$, and $T_w$ are the kinematic energy of the unimorph actuator, leading edge part (square bar on the top edge of the wing in Figure 1b), and wing, respectively. The total potential energy $U$ is also expressed by the sum of the elastic stored energy of the unimorph, $U_A$, and the pitching hinge, $U_w$, as

$$U=U_A+U_w.$$

(2)

By using $T$ and $U$ and the aerodynamic moment along the $\varphi$ and $\psi$ directions, ${\tau }_{\varphi }$ and ${\tau }_{\psi }$, respectively, the motion equation of the system can be written as

$$\begin{array}{c}\frac{\mathrm{d}}{\mathrm{d}t}\frac{\partial \Lambda }{\partial \dot{\varphi }}-\frac{\partial \Lambda }{\partial \varphi }={\tau }_{\varphi },\\ \frac{\mathrm{d}}{\mathrm{d}t}\frac{\partial \Lambda }{\partial \dot{\psi }}-\frac{\partial \Lambda }{\partial \psi }={\tau }_{\psi },\\ \Lambda =T-U,\end{array}\}$$

(3)

where $\Lambda$ is the Lagrangian of the system.

2.2. Mechanical Model

Here, we describe the derivation of $T$ and $U$ from the material and dimensional parameters of the system. Figure 2 shows the dimensions and coordinates of the system (details of the wing geometry are explained in a later section). The unimorph is designed as a trapezoid. It can reduce the stress concentration at the fixed end against inertial load and the aerodynamic force of the wing, and it can also reduce the risk of brittle fracture [24]. The cross point of the extended line of the unimorph side edges is designed to coincide with the center of the wing ($x=L_A+R/2$, where $R$ and $L_A$ are the lengths of the wing and unimorph, respectively). We assumed that the rotational axis of the stroke motion is placed at $x=L_A/2$. This assumption can be valid when the inertial and aerial forces of the wing are considered as concentrated forces at $x=L_A+R/2$. Of course, this is a very rough simplification, and it is not always suitable. Thus, we check the position of the rotational center of the wing experimentally; the result is described in the Result section. By using this assumption, the curvature of the unimorph $\kappa$ is simply expressed by $\varphi$ as

$$\kappa =\frac{\varphi }{L_A}.$$

(4)

The deflection of the unimorph $\delta$ can be given by

$$\delta =\frac{\kappa }{2}{x}^2.$$

(5)

Thus, $T_A$ is calculated from the following integral with respect to x:

$$T_A={{\displaystyle \int }}_0^{L_A}\frac{1}{2}\left({\displaystyle \sum }_i{\rho }_i{t}_i\right)w\left(x\right){\dot{\delta }}^2\mathrm{d}x=\frac{1}{240}\left({\displaystyle \sum }_i{\rho }_i{t}_i\right)\left(w_r+5w_t\right)L_A^3{\dot{\varphi }}^2,$$

(6)

where ${\rho }_i$ and $t_i$ are the density and thickness of the i-th layer of the unimorph, respectively, as shown in Figure 2. The bending moment caused by the piezoelectric effect is given by

$$M_p=E_0{d}_{pz}w\left(x\right)V\left(z_0-z_c\right),$$

(7)

where $E_0$ and $d_{pz}$ are Young’s modulus and the piezoelectric coefficient of the piezoelectric plate along the x-axis, respectively [25]. $V$ is the applied voltage to the piezoelectric material. $z_i$ are the z-positions of the layers of the unimorph, and they are expressed as $z_i={{\displaystyle \sum }}_j{t}_j-t_i/2$. $z_c$ is the neutral axis about the bending of the unimorph, and it is expressed as

$$z_c={\displaystyle \sum }_i{E}_i{t}_i{z}_i/{\displaystyle \sum }_i{E}_i{t}_i.$$

(8)

Assuming the unimorph as an Euler-Bernoulli beam [26], the curvature generated by the piezoelectric actuation is given by

$${\kappa }_{pz}=\frac{M_p}{S_{A}w\left(x\right)}=\frac{E_0{d}_{pz}V}{S_A}\left(z_0-z_c\right),$$

(9)

where $S_A$ is the bending stiffness of the unimorph per unit width, and it is given by

$$S_A={\displaystyle \sum }_i{E}_i\left(\frac{t_i^3}{12}+t_i{\left(z_i-z_c\right)}^2\right).$$

(10)

By using $\kappa$ and ${\kappa }_{pz}$, $U_A$ is calculated as

$$U_A={{\displaystyle \int }}_0^{L_A}\frac{1}{2}{S}_{A}w\left(x\right){\left(\kappa -{\kappa }_{pz}\right)}^2\mathrm{d}x=\frac{1}{2}{S}_A\frac{w_r+w_t}{L_A}{\varphi }^2-\frac{E_0{d}_{pz}V}{2}\left(z_0-z_c\right)\left(w_r+w_t\right).$$

(11)

Next, we describe the kinematic and potential energy of the wing. Figure 3 shows the coordinates and dimensional parameters. We newly define a coordinate $r$ as the distance from the wing root along the x-axis. $T_B$ is obtained as

$$T_B=\frac{1}{2}{J}_B{\dot{\varphi }}^2,$$

(12)

where $J_B$ is the inertia of the leading edge part around the stroke rotation axis. $T_w$ is calculated as

$$T_w=\frac{1}{2}{M}_w{\left|\mathit{v}\right|}^2+\frac{1}{2}{\mathit{\omega }}^T{\mathit{J}}_{\mathit{w}}\mathit{\omega },$$

(13)

where $M_w$ and ${\mathit{J}}_{\mathit{w}}$ are the mass and inertia tensor of the wing, respectively. The origin of ${\mathit{J}}_{\mathit{w}}$ is at the center of mass (COM) of the wing. $\mathit{v}$ and $\mathit{\omega }$ are the velocity of the COM and angular velocity of the wing, respectively. Upon defining these variables based on the wing-fixed x′y″z″ frame, $\mathit{v}$ and $\mathit{\omega }$ can be expressed by $\varphi$ and $\psi$ as

$$\mathit{v}=\left[\begin{array}{c}{v}_{x^{\prime }}\\ v_{y^{″}}\\ v_{z^{″}}\end{array}\right]=\left[\begin{array}{c}0\\ y_g\dot{\varphi }\sin \psi \\ \left(r_g+x_r\right)\dot{\varphi }+y_g\dot{\psi }\cos \psi \end{array}\right],$$

(14)

and

$$\mathit{\omega }=\left[\begin{array}{c}{\omega }_{x^{\prime }}\\ {\omega }_{y^{″}}\\ {\omega }_{z^{″}}\end{array}\right]=\left[\begin{array}{c}\dot{\psi }\\ \dot{\varphi }\cos \psi \\ \dot{\varphi }\sin \psi \end{array}\right],$$

(15)

where $y_g$ and $r_g$ are the $y$- and $r$-position of the COM, respectively, and $x_r$ is the distance between the wing root and the center axis of the stroke rotation. Upon assuming that the wing is a flat thin plate, the inertia terms with respect to x′z″ and y″z″ can be considered to be zero, and ${\mathit{J}}_{\mathit{w}}$ is expressed as

$${\mathit{J}}_{\mathit{w}}=\left[\begin{array}{ccc}{J}_{wx}& -J_{wxy}& 0\\ -J_{wxy}& J_{wy}& 0\\ 0& 0& J_{wx}+J_{wy}\end{array}\right]$$

(16)

By substituting Equations (14)–(16) into Equation (13), $T_w$ is rewritten as

$$\begin{gathered} T_w=\frac{1}{2}{M}_w\left[{\left(r_g+x_r\right)}^2{\dot{\varphi }}^2+y_g^2{\dot{\psi }}^2+2\left(r_g+x_r\right)y_g\dot{\varphi }\dot{\psi }\cos \psi \right]+\frac{1}{2}{J}_{wx}\left[{\dot{\psi }}^2+{\left(\dot{\varphi }\sin \psi \right)}^2\right] \\ +\frac{1}{2}{J}_{wy}{\dot{\varphi }}^2+J_{wxy}\dot{\varphi }\dot{\psi }\cos \psi . \end{gathered}$$

(17)

$k_w$ is defined as the spring constant of the pitching hinge, and $U_w$ is simply obtained as

$$U_w=\frac{1}{2}{k}_w{\psi }^2.$$

(18)

2.3. Aerodynamics on Wing

As stated above, we model the aerodynamic force as the sum of the translational force, force due to the added-mass effect, and rotational damping. The translational forces acting on a blade element of the wing (as shown in Figure 3) are given by [9]

$$\begin{array}{c}\mathrm{d}{F}_{\mathrm{tr},\mathrm{Drag}}=-\frac{1}{2}{\rho }_{\mathrm{Air}}{C}_D\left(\alpha \right)c\left(r\right)u_e^2 \mathrm{sign}\left(u_e\right) \mathrm{d}r,\\ \mathrm{d}{F}_{\mathrm{tr},\mathrm{Lift}}=\frac{1}{2}{\rho }_{\mathrm{Air}}{C}_L\left(\alpha \right)c\left(r\right)u_e^2 \mathrm{sign}\left(u_e\right) \mathrm{d}r,\end{array}\}$$

(19)

where $\mathrm{d}{F}_{\mathrm{tr},\mathrm{Drag}}$ and $\mathrm{d}{F}_{\mathrm{tr},\mathrm{Lift}}$ are the drag force (along z′-axis) and lift force (along y-axis), respectively. ${\rho }_{\mathrm{Air}}$ is the density of air. $C_D\left(\alpha \right)$ and $C_L\left(\alpha \right)$ are coefficients of the translational forces for $\mathrm{d}{F}_{\mathrm{tr},\mathrm{Drag}}$ and $\mathrm{d}{F}_{\mathrm{tr},\mathrm{Lift}}$, respectively. $\alpha$ is the angle of attack of the wing ($\alpha =90°-\psi$). $c\left(r\right)$ is the chord length of the wing. $u_e$ is the velocity at the center of a blade element in the z’-direction, and it is expressed as $u_e=\left(r+x_r\right)\dot{\varphi }-\frac{1}{2}c\left(r\right)\dot{\psi }\cos \psi$. The force component normal to the wing plane $\mathrm{d}{F}_{\mathrm{tr},\mathrm{N}}$ is also given by [9]

$$\begin{array}{c}\mathrm{d}{F}_{\mathrm{tr},\mathrm{N}}=-\frac{1}{2}{\rho }_{\mathrm{Air}}{C}_N\left(\alpha \right)c\left(r\right)u_e^2 \mathrm{sign}\left(u_e\right) \mathrm{d}r,\\ C_N\left(\alpha \right)=C_L\left(\alpha \right)\cos \alpha +C_D\left(\alpha \right)\sin \alpha .\end{array}\}$$

(20)

By integrating $\mathrm{d}{F}_{\mathrm{tr},\mathrm{Lift}}$ with respect to $r$, the total translational lift force $F_{\mathrm{tr},\mathrm{Lift}}$ can be calculated as

$$\begin{gathered} F_{\mathrm{tr},\mathrm{Lift}}={{\displaystyle \int }}_0^R\mathrm{d}{F}_{\mathrm{tr},\mathrm{Lift}} \\ =\frac{1}{2}{\rho }_{\mathrm{Air}}{C}_L\left(\alpha \right)A_w\left[I_{w2,1}{R}^2{\dot{\varphi }}^2-I_{w1,2}R\overline{c}\dot{\varphi }\dot{\psi }\cos \psi +\frac{1}{4}{I}_{w0,3}{\overline{c}}^2{\left(\dot{\psi }\cos \psi \right)}^2\right]\mathrm{sign}\left(\overline{u_e}\right), \end{gathered}$$

(21)

where $\overline{c}$ is the mean chord length ($\overline{c}=A_w/R$, $A_w$ is the area of the wing). $I_{wi,j}$ is a nondimensional geometrical parameter of the wing, and it is defined as

$$I_{wi,j}={{\displaystyle \int }}_0^1\widehat{c}{\left(\widehat{r}\right)}^j {\left(\widehat{r}+\widehat{x_r}\right)}^i\mathrm{d}\widehat{r},$$

(22)

where $\widehat{r}=r/R$, $\widehat{x_r}=x_r/R$, and $\widehat{c}=c/\overline{c}$. In the integration in Equation (21), we approximated $\mathrm{sign}\left(u_e\right)$ as being uniform for the entire wing; we replaced $\mathrm{sign}\left(u_e\right)$ with $\mathrm{sign}\left(\overline{u_e}\right)$, where $\overline{u_e}$ is the value of $u_e$ at $\mathrm{r}=R/2$. This approximation can become invalid when the wing spins around itself; however, we concluded that it is acceptable because such motion is rare in the intended flapping behavior.

The translational moment about the pitching (x’-axis) $M_{\mathrm{tr},x^{\prime }}$ is also derived by integration with respect to $r$ as

$$\begin{gathered} M_{\mathrm{tr},x^{\prime }}={{\displaystyle \int }}_0^R\mathrm{c}\left(r\right)d_{\mathrm{cp}}\left(\alpha \right) \mathrm{d}{F}_{\mathrm{tr},\mathrm{N}} \\ =-\frac{1}{2}{\rho }_{\mathrm{Air}}{C}_N\left(\alpha \right)A_w\overline{c} d_{cp}\left(\alpha \right)\left[I_{w2,2}{R}^2{\dot{\varphi }}^2-I_{w1,3}R\overline{c}\dot{\varphi }\dot{\psi }\cos \psi +\frac{1}{4}{I}_{w0,4}{\overline{c}}^2{\left(\dot{\psi }\cos \psi \right)}^2\right]\mathrm{sign}\left(\overline{u_e}\right), \end{gathered}$$

(23)

where $d_{cp}$ is the normalized center position of the air pressure. Dickson et al. reported that it can be estimated as [8]

$$d_{\mathrm{cp}}\left(\alpha \right)=\frac{0.82}{\pi }\left|\alpha \right|+0.05.$$

(24)

The translational moment about the stroke (y-axis) $M_{\mathrm{tr},y}$ is obtained in a similar way:

$$\begin{gathered} M_{\mathrm{tr},y}=-{{\displaystyle \int }}_0^R\left(r+x_r\right) \mathrm{d}{F}_{\mathrm{tr},\mathrm{Drag}}\mathrm{sign}\left(\overline{u_e}\right) \\ =-\frac{1}{2}{\rho }_{\mathrm{Air}}{C}_D\left(\alpha \right)A_{w}R\left[I_{w3,1}{R}^2{\dot{\varphi }}^2-I_{w2,2}R\overline{c}\dot{\varphi }\dot{\psi }\cos \psi +\frac{1}{4}{I}_{w1,3}{\overline{c}}^2{\left(\dot{\psi }\cos \psi \right)}^2\right]\mathrm{sign}\left(\overline{u_e}\right). \end{gathered}$$

(25)

The translational coefficients are calculated using the equations below [4,27]:

$$\begin{array}{c}{C}_L\left(\alpha \right)=\frac{1}{2}\pi \frac{R}{\overline{c}}{\left(1+\sqrt{{\left(\frac{1}{2}\frac{R}{\overline{c}}\right)}^2+1}\right)}^{-1}\sin 2\alpha ,\\ C_D\left(\alpha \right)=\frac{C_{D,\max }+C_{D,\min }}{2}-\frac{C_{D,\max }-C_{D,\min }}{2}\cos 2\alpha ,\end{array}\}$$

(26)

where we used $C_{D,\max }=5.0$ and $C_{D,\min }=0.4$ based on [4].

The force $\mathrm{d}{F}_{\mathrm{am}}$ due to the added-mass effect on a blade element is expressed as [5,6,7,8,9]

$$\mathrm{d}{F}_{\mathrm{am}}=-\pi {\rho }_{\mathrm{Air}}{\left(\frac{c\left(r\right)}{2}\right)}^2\dot{u_n} \mathrm{d}r,$$

(27)

where $u_n$ is the normal component of the wing velocity: $u_n=\left(r+x_r\right)\dot{\varphi }\cos \psi -\frac{1}{2}c\left(r\right)\dot{\psi }$. $\dot{u_n}$ is obtained by a time-derivative of $u_n$ as $\dot{u_n}=\left(r+x_r\right)\left(\ddot{\varphi }\cos \psi -\dot{\varphi }\dot{\psi }\sin \psi \right)-\frac{1}{2}c\left(r\right)\ddot{\psi }$. The total force $F_{\mathrm{am}}$ is obtained by integration as

$$\begin{gathered} F_{\mathrm{am}}={{\displaystyle \int }}_0^R\mathrm{d}{F}_{\mathrm{am}} \\ =-\frac{\pi }{4}{\rho }_{\mathrm{Air}}{\overline{c}}^2{R}^2{I}_{w1,2}\left(\ddot{\varphi }\cos \psi -\dot{\varphi }\dot{\psi }\sin \psi \right)+\frac{\pi }{8}{\rho }_{\mathrm{Air}}{\overline{c}}^{3}R I_{w0,3}\ddot{\psi }. \end{gathered}$$

(28)

$F_{\mathrm{am}}$ acts along the z″-axis (normal direction of the wing). The added-mass moment about the pitching x′-axis $M_{\mathrm{am},{\mathrm{x}}^{\prime }}$ and y″-axis $M_{\mathrm{am},{\mathrm{y}}^{″}}$ are respectively given by

$$\begin{gathered} M_{\mathrm{am},{\mathrm{x}}^{\prime }}=-{{\displaystyle \int }}_0^R\frac{1}{2}\mathrm{c}\left(r\right) \mathrm{d}{F}_{\mathrm{am}} \\ =\frac{\pi }{8}{\rho }_{\mathrm{Air}}{\overline{c}}^3{R}^2{I}_{w1,3}\left(\ddot{\varphi }\cos \psi -\dot{\varphi }\dot{\psi }\sin \psi \right)-\frac{\pi }{16}{\rho }_{\mathrm{Air}}{\overline{c}}^{4}R I_{w0,4}\ddot{\psi } \end{gathered}$$

(29)

and

$$\begin{gathered} M_{\mathrm{am},{\mathrm{y}}^{″}}={{\displaystyle \int }}_0^R\left(r+x_r\right) \mathrm{d}{F}_{\mathrm{am}} \\ =-\frac{\pi }{4}{\rho }_{\mathrm{Air}}{\overline{c}}^2{R}^3{I}_{w2,2}\left(\ddot{\varphi }\cos \psi -\dot{\varphi }\dot{\psi }\sin \psi \right)+\frac{\pi }{8}{\rho }_{\mathrm{Air}}{\overline{c}}^3{R}^2 I_{w1,3}\ddot{\psi }. \end{gathered}$$

(30)

The moment due to the rotational damping on a blade element is calculated as [4]

$$\mathrm{d}{M}_{\mathrm{rot}}=-\frac{1}{2}{\rho }_{\mathrm{Air}}{C}_{rd}{\dot{\psi }}^2{{\displaystyle \int }}_{-c\left(r\right)}^0\left|y^3\right|\mathrm{d}y \mathrm{sign} \dot{\psi },$$

(31)

where $C_{rd}$ is a constant coefficient. Different $C_{rd}$ values of 2.0–5.0 were used, as reported in several previous studies as well [4,28,29]. In this study, we used $C_{rd}$ = 3.5. Upon integrating $\mathrm{d}{M}_{\mathrm{rot}}$, the total moment $M_{\mathrm{rot}}$ can be expressed as

$$M_{\mathrm{rot}}=-\frac{1}{2}{\rho }_{\mathrm{Air}}{C}_{rd}R I_{rd} {\dot{\psi }}^2 \mathrm{sign} \dot{\psi },$$

(32)

where $I_{rd}$ is a constant that depends on the wing geometry as follows:

$$I_{rd}={{\displaystyle \int }}_0^1{{\displaystyle \int }}_{-c\left(r\right)}^0\left|y^3\right|\mathrm{d}y\mathrm{d}\widehat{r}.$$

(33)

${\tau }_{\varphi }$ and ${\tau }_{\psi }$ are respectively given by

$$\begin{array}{c}{\tau }_{\varphi }=M_{\mathrm{tr},y}+M_{\mathrm{am},{\mathrm{y}}^{″}}\cos \psi ,\\ {\tau }_{\psi }=M_{\mathrm{tr},x^{\prime }}+M_{\mathrm{am},{\mathrm{x}}^{\prime }}+M_{\mathrm{rot}}.\end{array}\}$$

(34)

The total lift force $F_{\mathrm{Lift}}$ is given by

$$F_{\mathrm{Lift}}=F_{\mathrm{tr},\mathrm{Lift}}+F_{\mathrm{am}}\sin \psi .$$

(35)

2.4. Electrical Model of Unimorph Actuator

The current flow into the actuator is given by a time-derivative of the sum of the charge due to the static capacitance, $q_{\mathrm{stat}}$, and the charge due to the dynamic piezoelectric effect, $q_{\mathrm{dyna}}$, as

$$j_A=\frac{\mathrm{d}}{\mathrm{d}t}\left(q_{\mathrm{stat}}+q_{\mathrm{dyna}}\right)$$

(36)

$q_{\mathrm{stat}}$ and $q_{\mathrm{dyna}}$ are respectively calculated by integrating charge densities along the unimorph longitudinal direction as [30]

$$q_{\mathrm{stat}}={{\displaystyle \int }}_0^{L_A}{ϵ}_{pz}\frac{w\left(x\right)}{t_0}V\mathrm{d}x={ϵ}_{pz}\frac{w_r+w_t}{2t_0}{L}_{A}V$$

(37)

and [30]

$$q_{\mathrm{dyna}}={{\displaystyle \int }}_0^{L_A}{d}_{pz}{E}_0\left(y_c-y_0\right)w\left(x\right)\kappa \mathrm{d}x=d_{pz}{E}_0\left(y_c-y_0\right)\frac{w_r+w_t}{2t_0}\varphi ,$$

(38)

where ${ϵ}_{pz}$ is the permittivity of the piezoelectric plate. By substituting Equations (37) and (38) into Equation (36), $j_A$ can be rewritten as

$$j_A={ϵ}_{pz}\frac{w_r+w_t}{2t_0}{L}_A\dot{V}+d_{pz}{E}_0\left(y_c-y_0\right)\frac{w_r+w_t}{2t_0}\dot{\varphi }$$

(39)

and the consumed power of the system is obtained as

$$P_A=j_{A}V.$$

(40)

By using the above equations (Equations (6), (11)–(13), (18), and (34)), the motion equation (Equation (3)) is expressed by known parameters and $V$. Therefore, by giving a certain initial condition and a prescribed $V$, the motion of the system can be calculated by integrating Equation (3). In this study, we numerically performed the integration using the built-in function for ordinary differential equations, NDSolve, in Mathematica (Wolfram). The total lift force $F_{\mathrm{Lift}}$ and consumed power $P_A$ are obtained by substituting the calculated motion, $\varphi$ and $\psi$, and $V$ into Equations (33) and (38).

3. Experiment

3.1. Design and Fabrication

To verify the analytical model, we measured the performance of prototypes of the flapping-wing actuation system and compared the result to the analytical calculation. Figure 4 shows the structure of the system. The system comprises two parts: a piezoelectric unimorph actuator and a wing. The wing comprises four parts (shown in Figure 4a): Leading-edge bar, flexible vein, reinforcing vein, and wing membrane. The wing membrane is a thin polyester film, and it is supported by the veins. The veins consist of a combination of titanium (Ti) reinforcing layer and polyimide (PI) flexible layer. The wing can rotate around the x-axis at the hinge structure, which is also made of the PI layer. A magnified drawing of the hinge structure is shown in Figure 4c. The hinge of the prototype comprises a torsional beam structure. Thus, $k_w$ can be calculated simply as the spring constant of the torsional beams [31]:

$$k_w=N\beta G_h\frac{w_h{t}_h^3}{L_h}, \beta \cong \frac{1}{3}-0.21\frac{t_h}{w_h}\left(1-\frac{1}{12}{\left(\frac{t_h}{w_h}\right)}^4\right),$$

(41)

where $N$ and $G_h$ are the number of torsional beams and the shear modulus of PI, respectively. $t_h$, $w_h$, and $L_h$ are the thickness, width, and length of the torsional beams, respectively. The unimorph actuator consists of a piezoelectric plate and a Ti shim plate bonded by epoxy adhesive (shown in Figure 4b). We used a Pb(In_1/2Nb_1/2)O₃-Pb(Mg_1/3Nb_2/3)O₃-PbTiO₃ (PIN-PMN-PT) single crystal as the piezoelectric material; it is suitable for large-displacement actuators because of its extremely high piezoelectric coefficient. PIN-PMN-PT plates were obtained from TRS Technologies (United States). The Ti shim plate has slits along the longitudinal direction to let excess adhesive flow out.

In this experiment, two prototypes with different size were measured for confirming the generality of the model. One has a 17.0-mm-long wing and a 10.5-mm-long unimorph and the other has a 32.4-mm-long wing and a 21.0-mm-long unimorph; they are respectively called “small type” and “large type” below. The other design parameters and material properties are summarized in Appendix A. Figure 5 shows photographs of the small and large types.

3.2. Measurement

Figure 6 shows the experimental setup for the wing motion, lift force, and power consumption. We captured the wing motion using an experimental setup with a high-speed camera. The average lift force and consumed power were measured using a precise electrical balance (HR-100A, A&D Co., Tokyo, Japan) and a power analyzer (PA1000, Tektronix, Beaverton, United States), respectively. To block the wind produced by the wing, an acrylic plate was placed between the wing and the electric balance.

Figure 7 shows a typical image of the actuated wing. We tracked the tip and root and the tip position of the leading edge, ${\mathit{p}}_{tip}$ and ${\mathit{p}}_{root}$, respectively, and a position on the lower edge of the wing, ${\mathit{p}}_{edge}$, based on the images to estimate $\varphi$ and $\psi$. $\varphi$ was simply calculated as the angle between ${\mathit{p}}_{tip}-{\mathit{p}}_{root}$ and the x-axis. $\psi$ was estimated by

$$\psi =\arcsin \frac{d}{c_{\mathrm{edge}}},$$

(42)

where $d$ and $c_{\mathrm{edge}}$ are the distance between the leading edge and ${\mathit{p}}_{edge}$ in the image and the chord length at the measured point, respectively. We applied unipolar sinusoidal driving signals with amplitude of $V_0$: $V\left(t\right)=V_0\sin \left({\omega }_{r}t\right)+V_0$. In this experiment, the maximum voltages for the small and large types were 30 and 50 V, respectively, to prevent fracture of the piezoelectric material. The frequency of the voltage signal ${\omega }_r$ was matched to the resonance frequency of the system.

4. Results and Discussion

4.1. Center Position of Stroke Motion

We estimated the rotational center position, ${\mathit{p}}_{center}$, of the wing to check whether approximation 1 is acceptable. The center positions were determined as the nearest position from lines through ${\mathit{p}}_{tip}$ and ${\mathit{p}}_{root}$. Specifically, we solved the following problem:

$$\underset{{\mathit{p}}_{center}}{\mathrm{minimize}}{\displaystyle \sum }_i{\left|\frac{\left({\mathit{p}}_{root}^{\mathit{i}}-{\mathit{p}}_{tip}^{\mathit{i}}\right)}{\left|{\mathit{p}}_{root}^{\mathit{i}}-{\mathit{p}}_{tip}^{\mathit{i}}\right|}\times \left({\mathit{p}}_{center}-{\mathit{p}}_{tip}^{\mathit{i}}\right)\right|}^2,$$

(43)

where ${\mathit{p}}_{root}^{\mathit{i}}$ and ${\mathit{p}}_{tip}^{\mathit{i}}$ are the positions in the i-th image frame. Figure 8 shows the tracking and estimation results; the blue solid lines indicate the line segments between ${\mathit{p}}_{root}^{\mathit{i}}$ and ${\mathit{p}}_{tip}^{\mathit{i}}$, and the red point indicates the estimated center position ${\mathit{p}}_{center}$. Both types show rotational motion centering on an almost fixed point. The positional errors relative to $L_A$ between ${\mathit{p}}_{center}$ and the assumed center position (center at the unimorph) were 4.1% for the small type and 8.4% for the large type. We concluded that these errors were small and acceptable to demonstrate the effectiveness of the model.

4.2. Resonance Frequency

The resonance frequencies were identified by the stroke amplitude response against the frequency-swept input signal in the experiment; we defined the resonant frequency as the frequency with the maximum stroke amplitude. In the calculation, we estimated the resonance frequencies from impulse system responses. The motion equation was solved with an initial condition of $\varphi =\dot{\varphi }=\psi =\dot{\psi }=0$ at $t=0$ and an impulse input of V. The resonant frequency ${\omega }_r$ can be derived by fitting $\varphi \left(t\right)$ to a function of ${\omega }_r$: $f\left(t\right)=A e^{-\lambda t}\sin \left({\omega }_{r}t+\theta \right)$, where $A$, $\lambda$, and $\theta$ are parameters that are also estimated in the fitting operation. This method has an advantage of shorter computational time comparing to calculating the response to the frequency-swept input. Figure 9 shows the calculated impulse responses and their fitted curves.

Figure 10 shows a comparison of the measured and calculated results. Blue and red bars show the measured and calculated data, respectively. The measure resonant frequencies of the small and large types were 210 and 105 Hz, respectively. The corresponding calculated resonant frequencies were 209 and 109 Hz. The small type showed roughly two times higher resonance frequency than the large type. The calculations were in very good agreement with the experiments; the errors of the small and large types were −0.7% and 3.8%, respectively. This result indicates that the kinematic part (inertias and springs) of the system has been modelled accurately.

4.3. Resonant Driving Characteristics

In this experiment, we measured the resonant driving performance of the small- and large-type actuators with $V_0$ of 15–30 and 20–50 V, respectively. Figure 11 and Figure 12 show the motion of the small and large types, respectively. Given the simplicity of the analytical model, we found that the behaviors of the calculated $\varphi$ and $\psi$ showed good agreement with the measured data. Figure 13 and Figure 14 show the amplitudes of $\varphi$ and $\psi$. The average absolute errors of $\varphi$ and $\psi$ amplitude were 20% and 8.1%, respectively. Figure 15 and Figure 16 show these tendencies of $\varphi$ with the lift force and power consumption, respectively. The average absolute errors of the lift force and consumed power between the calculation and the experiment were 15% and 8.7%, respectively. These results suggest that the derived model can well simulate the performances of the flapping-wing actuation system. Figure 17 shows the lift-to-power efficiency of the actuation system; it is defined as the ratio of the lift force to the power consumption. Upon increasing the applied voltage, the efficiency of both types of actuators decreased. The efficiencies of the small and large types at their maximum applied voltages were 10.8 and 7.31 gf/W, respectively.

There are two possible causes of lift force error: Wing motion ($\varphi$, $\psi$) and analytical model of lift force. As shown in Figure 13 and Figure 14, the predicted $\varphi$ and $\psi$ have nonnegligible error. In particular, the trend of the error in $\varphi$ showed a similarity to the lift force error. For example, for the small-type actuator, lower $V_0$ resulted in underestimated $\varphi$ (calculated $\varphi$ was smaller than the measured data) and larger $V_0$, in overestimated $\varphi$. On the other hand, for the large-type actuator, the calculated $\varphi$ was slightly smaller than the measured one for all tested $V_0$ values. The lift forces showed the same trends (Figure 15 and Figure 16). Thus, the error of the lift force was considered to be dominantly caused by the estimation error in the stroke motion. To verify this result, we performed additional calculations in which the lift force was estimated from the measured trajectory of $\varphi$ and $\psi$ using Equation (35). The results are shown in Figure 18 and Figure 19; in these figures, we compare the lift forces obtained using the experiment, calculation with the proposed aeromechanic model, and calculation with the prescribed measured trajectories. For both small- and large-type actuators, the lift force obtained from the measured trajectories resulted in smaller errors than that obtained from the fully analytical calculation. This result confirms that the prediction error of the wing motion is the main error source in the lift force estimation. We expect that the accuracy of the wing motion prediction is limited by the imprecision of the aerodynamic drag force and that improving the drag force model is key in the performance analysis of flapping-wing actuators.

5. Conclusions

We modeled a flapping-wing actuation system as a 2-DOF kinematic system with the semi-empirical quasi-steady aerodynamic forces and the electromechanical effect of piezoelectricity. We fabricated actuators of two different scales and measured their performances. The experimental results were in good agreement with the analytical calculations for both types of actuators, and the errors in the evaluated characteristics were less than 30%. We concluded that the analytical model effectively simulates the actual prototypes.