Distributed flow sensing for closed-loop speed control of a flexible fish robot

The fish robot modeled as a Joukowski foil takes the shape of the output of the Joukowski transformation of a circle. The Joukowski transformation, which is essentially a conformal mapping, is expressed as [14]

$$\begin{eqnarray}&&{z}^{0}=\xi +\displaystyle \frac{{a}^{2}}{\xi },\end{eqnarray} \tag{ 1 }$$

where the set of points ξ is prescribed to be a circle with radius R centered at ξ₀ in the complex ξ-plane. The x coordinate of the intersection of the circle and the x^ξ-axis is the Joukowski transformation parameter a, which is approximately one quarter of the chord length l of the foil. The image of the mapping in the z⁰-plane defines the boundary of the fish robot (figure 1). The origin of the z⁰-plane is O⁰. The x⁰-axis points along the chord line from the leading edge to the trailing edge.

Zoom In Zoom Out Reset image size

Define the body-fixed reference frame (the z-plane) by translating the z⁰-plane from point O⁰ to the point O, about which the fish robot flaps or rotates. The center of rotation for the flapping motion is chosen to be the one-quarter point along the camber line, as measured from the leading edge. Let z⁰_O be the coordinate of point O in the z⁰-plane. The transformation from the ξ-plane to the z-plane is

$$\begin{eqnarray}&&z=\xi +\displaystyle \frac{{a}^{2}}{\xi }-{z}_{O}^{0}.\end{eqnarray} \tag{ 2 }$$

A two-dimensional fluid with velocity U_f flows past the foil-shaped fish robot. The inertial frame I is defined so that the incoming flow velocity U_f is aligned with the x^I-axis. The flow velocity relative to the body is denoted U. The angle between the x-axis and the direction of the relative velocity U is the angle of attack α, with the nose pitching up chosen to be the positive direction (figure 1).

Assume that the circle radius R and the area of the Joukowski foil are constant during the flapping motion, even when the shape changes. (This assumption is based on the incompressibility of the foil material.) The area of the Joukowski foil [33] is

$$\begin{eqnarray}&&S=\pi {R}^{2}\left(1-\displaystyle \frac{{a}^{4}}{{({R}^{2}-{| {\xi }_{0}| }^{2})}^{2}}\right).\end{eqnarray} \tag{ 3 }$$

Parameter k describes the instantaneous shape of the fish robot, i.e.

$$\begin{eqnarray}&&k=\displaystyle \frac{{a}^{2}}{{R}^{2}-{| {\xi }_{0}| }^{2}}.\end{eqnarray} \tag{ 4 }$$

Along with the constant R, the shape of the fish robot is determined by the placement of the circle center ${\xi }_{0}={x}_{{\xi }_{0}}+{\rm{i}}{y}_{{\xi }_{0}},$ as described next.

In the foil plane, there are two shape parameters for the Joukowski foil: the camber ratio H and the thickness ratio T [14]. Intuitively, the camber ratio describes how much the foil bends. Under the assumption that $| {\xi }_{0}| \ll a,$ which is normally true for a moderately cambered Joukowski foil, H and T are linearly dependent on the vertical and horizontal displacement of the center of the circle, ${x}_{{\xi }_{0}}$ and ${y}_{{\xi }_{0}},$ respectively, i.e.

$$\begin{eqnarray}&&H=\displaystyle \frac{{y}_{{\xi }_{0}}}{2a},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&T=-\displaystyle \frac{3\sqrt{3}}{4}\displaystyle \frac{{x}_{{\xi }_{0}}}{a}.\end{eqnarray} \tag{ 6 }$$

Assume that the flexible fish robot changes its shape while maintaining a cambered Joukowski foil profile, meaning that each state of the deformation corresponds to a Joukowski foil shape with a varying camber ratio H. The camber ratio captures the degree of bending during the deformation. Meanwhile, given a camber ratio H that corresponds to the swimming state of the fish robot, the shape parameters in the ξ-plane, ${y}_{{\xi }_{0}},{x}_{{\xi }_{0}},$ and a, are functions of R, H and k.

According to potential-flow theory [14, 15], the flow in the complex circle plane (ξ-plane) generates the flow in the corresponding foil plane (z-plane), according to a conformal map. In an inviscid, incompressible, and irrotational fluid, the quasi-steady complex potential of the flow in the ξ-plane is a function of the relative flow speed U, the angle of attack α, the radius R, and the center ξ₀. The complex potential [14]

$$\begin{eqnarray}W(\xi ) & = & U(\xi -{\xi }_{0}){{\rm{e}}}^{-{\rm{i}}\alpha }+U\displaystyle \frac{{R}^{2}}{\xi -{\xi }_{0}}{{\rm{e}}}^{{\rm{i}}\alpha }\\ & & +{\rm{i}}\displaystyle \frac{{\rm{\Gamma }}}{2\pi }\mathrm{ln}(\xi -{\xi }_{0}),\end{eqnarray} \tag{ 7 }$$

represents the sum of three elementary flow fields: a uniform flow, a doublet, and a point vortex located at the center of the circle. The vortex circulation Γ is evaluated by enforcing the Kutta condition [14], which requires the trailing edge to be a stagnation point. The vortex circulation is [14]

$$\begin{eqnarray}&&{\rm{\Gamma }}=4\pi {RU}\mathrm{sin}(\alpha +\beta ).\end{eqnarray} \tag{ 8 }$$

Here $\beta =\mathrm{arcsin}({x}_{{\xi }_{0}}/R)\approx 2H$ is the phase angle of the center point ξ₀. Under the assumption that $| {\xi }_{0}| \ll a,\beta$ is approximated by 2H, which leads to

$$\begin{eqnarray}&&{\rm{\Gamma }}=4\pi {RU}\mathrm{sin}(\alpha +2H).\end{eqnarray} \tag{ 9 }$$

The conjugate flow velocity $\bar{f(z)}=u-{\rm{i}}v$ in the z-plane (the overline notation $\bar{\cdot }$ denotes the conjugate operator) is calculated using the complex potential in the ξ-plane and the Joukowski transform function, which yields

$$\begin{eqnarray}\bar{f(z)} & = & \displaystyle \frac{\partial W}{\partial \xi }{\left(\displaystyle \frac{\partial z}{\partial \xi }\right)}^{-1}\\ & = & \left(U{{\rm{e}}}^{-{\rm{i}}\alpha }-U\displaystyle \frac{{R}^{2}}{{(\xi -{\xi }_{0})}^{2}}{{\rm{e}}}^{{\rm{i}}\alpha }\right.\\ & & +\left.{\rm{i}}\displaystyle \frac{2{UR}\mathrm{sin}(\alpha +\beta )}{\xi -{\xi }_{0}}\right){\left(1-\displaystyle \frac{{a}^{2}}{{\xi }^{2}}\right)}^{-1},\end{eqnarray} \tag{ 10 }$$

where ξ is obtained using the inverse Joukowski transform

$$\begin{eqnarray}&&\xi =\xi (z)=0.5(z+{z}_{O}^{0})\\ &&\qquad \pm \sqrt{0.25{(z+{z}_{O}^{0})}^{2}-{a}^{2}}.\end{eqnarray} \tag{ 11 }$$

(The root outside the ξ-plane circle is used.) With the quasi-steady potential-flow model and the shape-parameter relationship in the Joukowski transformation, the flow field around the fish robot foil can be calculated given any parameter set (U, α, H).

The quasi-steady potential-flow model (10) does not describe the unsteady or transient effects caused by the flapping motion of a flexible fish robot. This subsection presents a second flow model used for simulation, which features discrete-time vortex shedding. In this model, a vortex is shed into the flow from the trailing edge of the foil at every discrete time step. The shed vortices convect with the flow according to the local fluid velocity. The net circulation around the robot body contributes to drag, which is not predicted by the quasi-steady potential-flow theory. The accuracy of the vortex-shedding method has been investigated and validated [33–35]. However, the system complexity increases with time and grows too fast for real-time use. We instead use the vortex-shedding method for simulating the flow field to test the effectiveness of the quasi-steady potential-flow model for flow-parameter estimation.

Let Ω be the angular velocity of the fish robot with counter-clockwise rotation about the pivot point O chosen to be the positive direction ($\dot{\alpha }=-{\rm{\Omega }}$). In the vortex-shedding model, the complex flow potential with respect to the ξ-plane is

$$\begin{eqnarray}W & = & U(\xi -{\xi }_{0}){{\rm{e}}}^{-{\rm{i}}\alpha }+U\displaystyle \frac{{R}^{2}}{\xi -{\xi }_{0}}{{\rm{e}}}^{{\rm{i}}\alpha }+{\rm{\Omega }}{W}_{{\rm{\Omega }}}\\ & & +{\rm{i}}\displaystyle \frac{{{\rm{\Gamma }}}_{0}}{2\pi }\mathrm{ln}\left(\displaystyle \frac{\xi -{\xi }_{0}}{R}\right)\\ & & +\displaystyle \sum _{k=1}^{n}{\rm{i}}\displaystyle \frac{{{\rm{\Gamma }}}_{k}}{2\pi }\left(\mathrm{ln}(\xi -{\xi }_{k})\right.\\ & & -\left.\mathrm{ln}\left(\xi -{\xi }_{0}-\displaystyle \frac{{R}^{2}}{\bar{{\xi }_{k}-{\xi }_{0}}}\right)\right).\end{eqnarray} \tag{ 12 }$$

The corresponding (unit) complex potential is denoted by W_Ω. Γ₀ represents the vortex circulation at the center of the circle and Γ_k represents the circulation of the kth vortex located at position ξ_k. The complex potential is

$$\begin{eqnarray}{W}_{{\rm{\Omega }}} & = & -\displaystyle \frac{{\rm{i}}}{2}\left({R}^{2}+2{\xi }_{0}\displaystyle \frac{{R}^{2}}{\xi -{\xi }_{0}}+{| {\xi }_{0}| }^{2}\right.\\ & & +2{a}^{2}\displaystyle \frac{\frac{{R}^{2}}{\xi -{\xi }_{0}}+\bar{{\xi }_{0}}}{\xi }+2{z}_{O}^{0}\displaystyle \frac{{R}^{2}}{\xi }+2\bar{{z}_{O}^{0}}({\xi }_{0}+\displaystyle \frac{{a}^{2}}{\xi })\\ & & +\left.{z}_{O}^{0}\bar{{z}_{O}^{0}}+\displaystyle \frac{{a}^{4}(\xi -2{\xi }_{0})}{\xi ({R}^{2}-{| {\xi }_{0}| }^{2})}\right).\end{eqnarray} \tag{ 13 }$$

The complex velocity in the z-plane using the vortex-shedding method is

$$\begin{eqnarray}&&\bar{f(z)}=\displaystyle \frac{\partial W}{\partial \xi }{\left(\displaystyle \frac{\partial z}{\partial \xi }\right)}^{-1}.\end{eqnarray} \tag{ 14 }$$

Equation (14) describes the entire flow field except at the vortex locations, which are singular points and thus the flow is undefined. Using Routh's rule [36, 37], the kth vortex conjugate velocity $\bar{f({z}_{k})}={u}_{{z}_{k}}-{\rm{i}}{v}_{{z}_{k}}$ is

$$\begin{eqnarray}\bar{f({z}_{k})} & = & \bar{\displaystyle \frac{{\rm{d}}{z}_{k}}{{\rm{d}}t}}=\displaystyle \frac{{\rm{d}}}{{\rm{d}}z}{\left.\left(W-{\rm{i}}\displaystyle \frac{{{\rm{\Gamma }}}_{k}}{2\pi }\mathrm{ln}(z-{z}_{k})\right)\right|}_{z={z}_{k}}\\ & = & {\left(\displaystyle \frac{{\rm{d}}z}{{\rm{d}}\xi }\right)}^{-1}\displaystyle \frac{{\rm{d}}}{{\rm{d}}\xi }\left(W(\xi )-{\rm{i}}\displaystyle \frac{{{\rm{\Gamma }}}_{k}}{2\pi }\mathrm{ln}(\xi -{\xi }_{k})\right)\\ & & -{\left.{\rm{i}}\displaystyle \frac{{{\rm{\Gamma }}}_{k}}{4\pi }\displaystyle \frac{{{\rm{d}}}^{2}z}{{\rm{d}}{\xi }^{2}}{\left(\displaystyle \frac{{\rm{d}}z}{{\rm{d}}\xi }\right)}^{-2}\right|}_{\xi ={\xi }_{k}}.\end{eqnarray} \tag{ 15 }$$

The circulation strength Γ_k of the vortex shed into the flow at the point ξ_k is evaluated by enforcing the Kutta condition at the trailing edge.

The position of the shed vortex z_k may be modeled using one of several existing approaches. This paper adopts the '1/3 arc-length' method introduced by Streitlien and Triantafyllou [38], which places the shed vortex at the one-third point of the arc tangent to the camber line at the trailing edge and passing through the previous vortex point: i.e.

$$\begin{eqnarray}{z}_{k} & = & ({z}_{k-1}-{z}_{{\rm{TE}}})\displaystyle \frac{1-{{\rm{e}}}^{{\rm{i}}\frac{\theta }{3}}}{1-{{\rm{e}}}^{{\rm{i}}\theta }}\\ & = & \displaystyle \frac{{z}_{k-1}-{z}_{{\rm{TE}}}}{1+{{\rm{e}}}^{{\rm{i}}\frac{\theta }{3}}+{{\rm{e}}}^{{\rm{i}}\frac{2\theta }{3}}}+{z}_{{\rm{TE}}}.\end{eqnarray} \tag{ 16 }$$

Here the trailing-edge coordinate is ${z}_{{\rm{TE}}}=2a-{z}_{O}^{0}$ and θ is the arc angle satisfying

$$\begin{eqnarray}&&{{\rm{e}}}^{{\rm{i}}\theta }={{\rm{e}}}^{-{\rm{i}}2(-2\beta )}\displaystyle \frac{{({z}_{k-1}-{z}_{{\rm{TE}}})}^{2}}{{\left|{z}_{k-1}-{z}_{{\rm{TE}}}\right|}^{2}},\end{eqnarray} \tag{ 17 }$$

where $\beta =\mathrm{arcsin}({y}_{{\xi }_{0}}/R)$ as before.

When the fish robot rotates about the point O, the z-reference frame also rotates with it, while the shed vortices move as predicted by Routh's rule in the inertial frame. Therefore, the vortex coordinates with respect to the z-plane are

$$\begin{eqnarray}&&{z}_{k}(N+1)=({z}_{k}(N)+f({z}_{k})T){{\rm{e}}}^{{\rm{i}}{\rm{\Delta }}\alpha },\end{eqnarray} \tag{ 18 }$$

where z_k(N) is the kth vortex coordinate in the z-plane at the Nth time step, T is the time step, and Δα is the change in the angle of attack from the Nth time step to the (N + 1)th time step.

Let ${\boldsymbol{\Pi }}={\bf{J}}\omega$ and ${\bf{P}}={\bf{Mv}}$ denote the total angular and linear momenta of the robot-fluid system, respectively, where ω is the angular velocity of the body-fixed frame with respect to the inertial frame expressed in the body-fixed frame, ${\bf{v}}$ is the corresponding translational velocity, ${\bf{J}}$ is the inertia matrix, and ${\bf{M}}$ is the mass matrix. Both ${\bf{J}}$ and ${\bf{M}}$ include the added-mass effect that describes the additional effect (force) resulting from water acting on the fish robot during acceleration or deceleration. Assume the off-diagonal terms in the inertial matrix are negligible and the added mass and added inertia do not vary significantly during the flapping motion. The dynamics of the fish robot are governed by Kirchhoff's equations [28, 39, 40], i.e.

$$\begin{eqnarray}&&\dot{{\boldsymbol{\Pi }}}={\boldsymbol{\Pi }}\times \omega +{\bf{P}}\times {\bf{v}}+{\bf{T}},\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&\dot{{\bf{P}}}={\bf{P}}\times \omega +{\bf{F}},\end{eqnarray} \tag{ 20 }$$

where ${\bf{T}}$ is the external moment vector and ${\bf{F}}$ is the external force vector.

For the planar motion of the fish robot, which is the focus of this paper, we have

$$\begin{eqnarray*}&&\omega ={[0,\;0,{\rm{\Omega }}]}^{T},\quad {\bf{v}}={[{v}_{1},{v}_{2},0]}^{T},\\ &&{\boldsymbol{\Pi }}={[0,\;0,J{\rm{\Omega }}]}^{T},\quad {\bf{P}}={[{m}_{1}{v}_{1},{m}_{2}{v}_{2},0]}^{T},\\ &&{\bf{T}}={[0,\;0,{T}_{{\rm{p}}}]}^{T},\quad {\rm{and}}\;{\rm{}}{\bf{F}}=\left[-\left({F}_{{\rm{t}}}-{F}_{{\rm{d}}}\mathrm{cos}\alpha \right.\right.\\ &&\quad +{\left.\left.{F}_{{\rm{l}}}\mathrm{sin}\alpha \right),\quad {F}_{{\rm{l}}}\mathrm{cos}\alpha +{F}_{{\rm{d}}}\mathrm{sin}\alpha ,0\right]}^{T}.\end{eqnarray*}$$

Here J is the sum of the inertia of the robot and the added inertia in the pitching direction, and m₁ and m₂ are the sum of the mass of the robot and the added mass in the surge and sway directions, respectively. Figure 2 illustrates the hydrodynamic pitching torque T_p; the thrust force F_t, generated by the flapping motion of the foil with the −x-axis direction as positive; the drag force F_d, in the opposite direction of the motion of the robot relative to the fluid; and the lift force F_l, perpendicular to the relative-flow direction.

Zoom In Zoom Out Reset image size

The dynamics for planar motion in first-order form are

$$\begin{eqnarray}&&J\dot{{\rm{\Omega }}}=({m}_{1}-{m}_{2}){v}_{1}{v}_{2}+{T}_{{\rm{p}}}\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{m}_{1}\dot{{v}_{1}}={m}_{2}{v}_{2}{\rm{\Omega }}-{F}_{{\rm{t}}}-{F}_{{\rm{d}}}\mathrm{cos}\alpha +{F}_{{\rm{l}}}\mathrm{sin}\alpha \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{m}_{2}\dot{{v}_{2}}=-{m}_{1}{v}_{1}{\rm{\Omega }}+{F}_{{\rm{l}}}\mathrm{cos}\alpha +{F}_{{\rm{d}}}\mathrm{sin}\alpha .\end{eqnarray} \tag{ 23 }$$

The hydrodynamic forces and moment are modeled following aerospace engineering conventions [41, 42], i.e.

$$\begin{eqnarray}&&{T}_{{\rm{p}}}={C}_{{\rm{p}}}(\alpha +2H){U}^{2}-{K}_{{\rm{p}}}{\rm{\Omega }}\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{F}_{{\rm{d}}}=({C}_{{\rm{d}}}^{0}+{C}_{{\rm{d}}}{(\alpha +2H)}^{2}){U}^{2}\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{F}_{{\rm{l}}}={C}_{{\rm{l}}}(\alpha +2H){U}^{2},\end{eqnarray} \tag{ 26 }$$

where ${C}_{{\rm{p}}},{C}_{{\rm{d}}}^{0},{C}_{{\rm{d}}},$ and C_l are hydrodynamic coefficients that can be identified using flow tunnel experiments, and K_p is the pitch damping coefficient.

In this paper, the angle of attack is the control input. Selecting a sinusoidal waveform for the control profile implies

$$\begin{eqnarray}&&{\alpha }_{\mathrm{ctrl}}=A\mathrm{sin}{\phi }_{\mathrm{ctrl}}=A\mathrm{sin}(2\pi {ft}),\end{eqnarray} \tag{ 27 }$$

where ϕ_ctrl is the phase angle of the angle of attack, and A and f represent the amplitude and frequency of the periodic control input α_ctrl, respectively. The thrust force generated by the periodic actuation (27) is approximated as [6]

$$\begin{eqnarray}&&{F}_{{\rm{t}}}={\bar{F}}_{{\rm{t}}}+({\hat{F}}_{{\rm{t}}}-{\bar{F}}_{{\rm{t}}})\mathrm{sin}(2{\phi }_{\mathrm{ctrl}}),\end{eqnarray} \tag{ 28 }$$

where ${\bar{F}}_{{\rm{t}}}$ and ${\hat{F}}_{{\rm{t}}}$ are the mean and maximum thrust force in one flapping period, respectively. The mean and maximum thrust force depend on the product of the amplitude and frequency of the flapping motion [6], i.e.

$$\begin{eqnarray}&&{\bar{F}}_{{\rm{t}}}={k}_{1}{({Af})}^{{k}_{2}},\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray}&&{\hat{F}}_{{\rm{t}}}={k}_{3}{({Af})}^{{k}_{4}}.\end{eqnarray} \tag{ 30 }$$

The parameters k₁, k₂, k₃, and k₄ are identified by force-sensing experiments (and typically k₂ and k₄ are approximately equal to 2 [6]).

Define V as the speed of the fish robot with respect to the inertial frame, with the −x^I-axis direction as positive. The inertial flow speed along the x^I-axis is U_f. Therefore, the relative speed of the flow with respect to the robot is U = U_f + V. We also have ${v}_{1}=-V\mathrm{cos}\alpha$ and ${v}_{2}=-V\mathrm{sin}\alpha .$

This paper focuses on one-dimensional swimming. The robot moves along the x^I-axis and the angle of attack is directly controlled, so the dynamics become

$$\begin{eqnarray}-{m}_{1}{m}_{2}\dot{U} & = & ({m}_{1}^{2}-{m}_{2}^{2})(U-{U}_{{\rm{f}}}){\rm{\Omega }}\mathrm{sin}\alpha \mathrm{cos}\alpha \\ & & -({F}_{{\rm{t}}}-{F}_{{\rm{d}}}\mathrm{cos}\alpha +{F}_{{\rm{l}}}\mathrm{sin}\alpha ){m}_{2}\mathrm{cos}\alpha \\ & & +({F}_{{\rm{l}}}\mathrm{cos}\alpha +{F}_{{\rm{d}}}\mathrm{sin}\alpha ){m}_{1}\mathrm{sin}\alpha .\end{eqnarray} \tag{ 31 }$$

Modeling camber dynamics is a challenging fluid-structure-interaction problem that may involve continuum mechanics and boundary-value partial differential equations. However, a tractable model captures the camber motion sufficiently well for the purpose of real-time control. The camber kinematics are a linear function of the second derivative of the angle of attack with respect to time, i.e.

$$\begin{eqnarray}&&\dot{H}=-{K}_{h}\ddot{\alpha }={K}_{h}\dot{{\rm{\Omega }}},\end{eqnarray} \tag{ 32 }$$

where K_h is the camber-dynamics coefficient.

The fish robot dynamics with sinusoidal actuation input is inspired by the flapping motion of swimming fish. The swimming speed U oscillates in a periodic waveform. For one-dimensional swimming with the angle of attack as the control input, we apply the averaging method [43] to capture the averaged dynamics of the motion, which is later adopted in the feedforward controller design described in section 4.

Rewrite the system equation (31) with periodic control input (27) in the state-space form

$$\begin{eqnarray}&&\dot{U}=f(t,U,u),\end{eqnarray} \tag{ 33 }$$

where u = [A, f]^T represents the parameterized control variables. The dynamic equation (33) is periodic with time period T = 1/f, i.e.

$$\begin{eqnarray}&&f(t,U,u)=f(t+T,U,u).\end{eqnarray} \tag{ 34 }$$

In order to understand the average effect of the oscillatory swimming motion, for example, to study the mean swimming speed, an averaging method is applied to obtain the average system [43]

$$\begin{eqnarray}&&\dot{U}={f}_{\mathrm{av}}(U,u),\end{eqnarray} \tag{ 35 }$$

where

$$\begin{eqnarray}&&{f}_{\mathrm{av}}(U,u)=\displaystyle \frac{1}{T}{\displaystyle \int }_{0}^{T}f(\tau ,U,u){\rm{d}}\tau .\end{eqnarray} \tag{ 36 }$$

To facilitate this calculation, assume that $\mathrm{sin}\alpha \approx \alpha$ and $\mathrm{cos}\alpha \approx 1,$ which is valid when α is sufficiently small. We also know that ${\int }_{0}^{T}\mathrm{sin}(4\pi f\tau ){\rm{d}}\tau =0$ and ${\int }_{0}^{T}\mathrm{cos}(4\pi f\tau ){\rm{d}}\tau =0.$ The average model (35) is expanded using (36) to obtain

$$\begin{eqnarray}&&\dot{U}=\displaystyle \frac{1}{{m}_{1}{m}_{2}}\left({m}_{2}{k}_{1}{({Af})}_{2}^{k}-g(A,f){U}^{2}\right),\end{eqnarray} \tag{ 37 }$$

where

$$\begin{eqnarray}g(A,f) & = & {m}_{2}\left({C}_{D}^{0}+{C}_{D}{A}^{2}/2+{C}_{L}{A}^{2}/2\right.\\ & & +\left.8{C}_{D}{\pi }^{2}{K}_{H}^{2}{({Af})}^{2}\right)\\ & & +{m}_{1}A\left({C}_{L}A/2+{C}_{D}^{0}A/2+3{C}_{D}{A}^{3}/8\right.\\ & & +\left.2{C}_{D}A{\pi }^{2}{K}_{H}^{2}{({Af})}^{2}\right).\end{eqnarray} \tag{ 38 }$$

The average model (37) requires the averaged system state $\frac{1}{T}{\int }_{0}^{T}U(\tau ){\rm{d}}\tau$ to change more slowly than the oscillating control input.

Silicone rubber (Ecoflex 00-30 from Smooth-On) is the flexible material chosen for the fish robot. This platinum-catalyzed type of silicone rubber is versatile and easy to use. Cured rubber is soft and strong with Shore 00-30 hardness, which rebounds to its original form without distortion after stretching. A mold of the fish robot was designed in SolidWorks, as shown in figures 7(a) and (b), and then the mold was manufactured using a high-precision 3D printer. The mixed silicone rubber compound is poured into the mold and then kept in a vacuum chamber during curing to avoid bubbles. For actuation of the angle of attack, a mini shaft from MakerBeam was inserted before the molding process at the one-quarter-point of the chord, behind the leading edge. For the pressure sensor ports, holes of radius 1 mm were placed on each side of the mold at a distance along the center line of 2.5, 7.5, and 10 cm behind the leading edge. We packaged six pressure sensors (Servoflo MS5401-BM) with appropriately sized tubing, and fixed the sensors to the center shaft. Each pressure sensor outputs analog voltage in proportion to the local pressure. The whole fish robot measures 20 cm long, 3.6 cm wide, and 12 cm tall as shown in figure 7(c).

Zoom In Zoom Out Reset image size

The experimental testbed for holding the fish robot and guiding the motion was designed to facilitate the flow-sensing and control experiment; a schematic is shown in figure 8. A double-rail linear guide is mounted to a 80/20 frame, which is placed over a Loligo flow tank. We used two linear motion air bearings (NEWWAY porous air bearing) to allow for low-friction horizontal motion. A gantry platform was mounted on those bearings to support the fish robot. Two load cell sensors (300 g measurement range) mounted vertically at equal distance from the center detect the force and torque in the surge and pitching directions, respectively. A dc servo (Savox 0235MG) sits between the two load cells. The servo arm is bolted to an L-shape MakerBeam rigidly attached to the fish to control the angle of attack. This setup was used for water-tunnel experiments to identify the hydrodynamic coefficients in table 1, and to investigate the thrust-force generation from the periodic flapping motion [6, 49]. When the air bearing is turned on, the fish robot moves freely along the guiding rail. The entire test platform is placed over a 185 L flow tank, which generates a uniform flow field within an enclosed test section that measures 25 × 25 × 87.5 cm, as shown in figure 9.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This subsection presents results from testing the recursive Bayesian estimation algorithm during one-dimensional swimming. The fish robot is still for the first 5 s, and then flaps for 10 s. The angle-of-attack actuation signal is sinusoidal with a frequency of 0.75 Hz and a constant amplitude of 10°. The uniform flow speed is 18.74 cm s⁻¹, which is calibrated using a OTT MF-Pro flow meter. The air bearing operates throughout the experiment to allow the fish to move forward and backward. The measurement data is acquired from the pressure sensors using DAQ 6225 from National Instruments. The data is then transmitted via USB to a laptop that runs the Bayesian filter for data assimilation in Matlab 2013b. The control commands for the angle of attack are sent via serial communication to an Arduino UNO that drives the servo. The system parameters identified before the experiments have the same values as listed in table 1. Figure 8 shows the data stream for the flow-sensing experiment.

Figure 10 presents the measurement data from all six pressure sensors. The data for each sensor has been pre-processed by subtracting the still-water pressure measurement taken before the experiment, in order to eliminate the influence of nonuniform and possibly time-varying bias from the pressure sensors. We find that the second and third sensors on each side have the peak pressure almost at the same time, whereas the first pressure sensor has a phase lag about 54° (0.2 s in time). Ideally, in a quasi-steady laminar flow chamber, all three sensors on the same side should have the same temporal behaviors. However, in our experimental setup, due to the size limitation of the flow tank, water is reflected from the wall of the tank and bounced back to different parts of the fish body at different times, which contributes to the phase difference. In addition, the first sensor rotates in the opposite direction compared to the other two sensors, which enhances the phase difference.

Zoom In Zoom Out Reset image size

Figure 11 shows the time evolution of the marginal probability densities. Observe that the estimated angle of attack α (figure 11(b)) has the same oscillation frequency as the actual angle of attack, with a slightly larger amplitude and a negative phase shift. The estimated camber ratio H (figure 11(c)) has a sinusoidal-like waveform, as predicted in simulation. The estimated relative speed U (figure 11(a)) oscillates around a constant value at twice the actuation frequency. The speed oscillation arises from the periodic thrust force generated by the flapping motion. The fish robot has the same inertial speed as the uniform flow and consequently exhibits station-holding behavior. (The gantry air tubing and electrical wiring restrict motion of the robot to within the test section.) The average estimated relative flow speed 18.74 cm s⁻¹ is equal to the incoming uniform flow speed, as observed in figure 11(a).

Zoom In Zoom Out Reset image size

Figure 12 shows the moving average of the estimated relative flow speed U, which is used as the observed true speed to be compared with the referenced value in the closed-loop control experiment. The moving-average window size is equal to one flapping time period. The shaded area shows the standard deviation of the estimate over the corresponding average time window. The results suggest that it takes only one flapping period for the averaged estimated flow-speed to converge, which permits fast closed-loop speed control performance. Note the flow-speed estimation is more accurate during fish flapping (5–15 s) than during fish resting (0–5 s), with significant difference in terms of estimation error (5% versus 60%). This observation shows that the fish robot senses unsteady flow better than steady flow at zero angle of attack; checking the observability Gramian of the robot's dynamic model with pressure measurements as the system output shows that the flow speed at zero angle of attack is unobservable [44, 46].

Zoom In Zoom Out Reset image size

This subsection describes the closed-loop speed-control experiment for a fish robot swimming freely in the surge direction. The objective is to match the speed between the robot locomotion and the incoming flow by controlling the flapping amplitude.

The fish robot rests in water at the beginning of the experiment for 5 s, and then conducts closed-loop speed control for 20 s. The flow tank speed starts at 20.6 cm s⁻¹ and then drops to 12.48 cm s⁻¹ at t = 16 s. The closed-loop control strategy combining feedforward and feedback control (see section 4) is implemented using distributed flow estimation to automatically regulate the flapping amplitude A. The actuation frequency is fixed at 0.75 Hz. The feedback control is calculated at the end of each flapping time period based on the average estimated flow speed during the previous interval.

Figure 13 shows the pressure-measurement data for the six sensors distributed equally on each side. The amplitude of each pressure oscillation changes significantly at t = 5 s and t = 16 s, which corresponds to the step change in the flow speed. At t = 5 s, the robot starts to flap, thus generating periodic pressure oscillations at each sensor location. At t = 16 s, the amplitude of pressure measurements drops significantly corresponding to the decrease in relative flow speed. A phase shift is observed between the first sensor on each side and the other two sensors on the same side as in the open-loop flow-sensing experiments.

Zoom In Zoom Out Reset image size

Figure 14 shows the time evolution of the marginal probability densities in the recursive flow estimation. The estimated flow speed oscillates around the desired reference speed and has a transient at t = 16 s when the flow speed decreases. The estimated angle of attack and camber ratio also reflect the step change. There is a noticable phase lag in angle-of-attack estimation as time progresses. Such an estimation error comes from the usage of a quasi-steady potential flow model in the recursive Bayesian filter, which ignores the unsteady effects. However, this simplification is necessary for the feasibility of real-time computation and control. In addition, it does not influence the control performance which depends on the average of the estimated angle of attack rather than the instantaneous value.

Zoom In Zoom Out Reset image size

Figure 15 shows the moving average of the estimated relative flow speed U with a time window size of one control/flapping period. The shaded area represents one standard deviation of U over the corresponding moving-average time window, illustrating the amplitude of the relative-flow speed oscillation. The average flow speed converges to the reference value. The transient from the step change in swimming speed lasts about three control periods, and the steady-state tracking error is 1 cm s⁻¹ (5%).

Zoom In Zoom Out Reset image size

Figure 16 shows the commanded angle of attack. The output of the feedforward controller is close to the total control output in the steady-state period, which indicates that the inverse mapping of the derived average model is a good choice for feedforward control. The tracking error is compensated by the feedback PI controller. The quick convergence and small tracking error further validate the control strategy and flow-sensing algorithm.

Zoom In Zoom Out Reset image size