Flow measurement and thrust estimation of a vibrating ionic polymer metal composite

Figure 4 shows the time trace of the tip displacement of the IPMC at the driving frequency $f=1.5\;\mathrm{Hz}$. Raw motion data are fitted with a sine-wave by utilizing the Matlab built-in function 'fit,' to demonstrate that the tip displacement is well described by a harmonic signal at the input frequency. In figure 5(a), we display the peak tip displacement ${\delta }_{p}$ as a function of the input frequency. The peak tip displacement decreases as the input frequency increases, likely due to the hydrodynamic damping from the encompassing fluid [54]. The hydrodynamic damping tends to mask the underwater resonance of the IPMC, which can be seen from the peak tip velocity in figure 5(b). Therein, we note that the fundamental resonance frequency is located around approximately $2.5\;\mathrm{Hz}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To interpret the experimental findings, we model the IPMC vibration by combining our recent physics-based modeling framework for IPMC actuation [7] with the classical Morison formula [42, 43] that describes fluid-structure interactions for slender bodies. Within this framework, the IPMC vibration is governed by the following equation:

$$\begin{eqnarray}&&K\displaystyle \frac{{\partial }^{4}w(x,t)}{\partial {x}^{4}}+\rho \displaystyle \frac{{\partial }^{2}w(x,t)}{\partial {t}^{2}}=H(x,t)+s(x,t),\end{eqnarray} \tag{ 1 }$$

where $w(x,t)$ is the deflection ($w(L,t)=\delta (t)$), $H(x,t)$ is the hydrodynamic forcing from the encompassing fluid, and $s(x,t)$ is the structural damping [55].

The boundary conditions are as follows [42, 56]:

$$\begin{eqnarray}&&w(0,t)=0,\end{eqnarray} \tag{ 2a }$$

$$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial w(x,t)}{\partial x}\right|}_{x=0}=0,\end{eqnarray} \tag{ 2b }$$

$$\begin{eqnarray}&&{\left.\displaystyle \frac{{\partial }^{3}w(x,t)}{\partial {x}^{3}}\right|}_{x=L}=0,\end{eqnarray} \tag{ 2c }$$

$$\begin{eqnarray}&&K{\left.\displaystyle \frac{{\partial }^{2}w(x,t)}{\partial {x}^{2}}\right|}_{x=L}=M(t),\end{eqnarray} \tag{ 2d }$$

where M(t) is the bending moment generated by counterion migration within the IPMC. Such a bending moment is associated with ion mixing and electric polarization within the ionomer core and the electrodes [7]. We hypothesize that the time scale of the IPMC deformation is considerably slower than the chemoelectric response so that the bending moment M(t) can be described from the quasi-static solution presented in [7], which posits a simple functional dependence of M(t) on the square of the applied voltage, without a significant dependence on the mechanical deformation. As a result, we simply treat the bending moment as a known function of time of the form $M(t)={M}_{p}\mathrm{sin}(\omega t)$.

Following [42, 43], we write the hydrodynamic loading $H(x,t)$ as

$$\begin{eqnarray}H(x,t) & = & -\displaystyle \frac{\pi }{4}{\rho }_{{\rm{f}}}{b}_{0}^{2}{C}_{{\rm{A}}}\displaystyle \frac{{\partial }^{2}w(x,t)}{\partial {t}^{2}}\\ & & -\displaystyle \frac{\pi }{2}{\rho }_{{\rm{f}}}{b}_{0}{C}_{{\rm{D}}}\displaystyle \frac{\partial w(x,t)}{\partial t}\left|\displaystyle \frac{\partial w(x,t)}{\partial t}\right|,\end{eqnarray} \tag{ 3 }$$

where ${C}_{{\rm{A}}}$ is the added mass coefficient and ${C}_{{\rm{D}}}$ is the hydrodynamic damping coefficient.

We approximate the structural damping with a simple hysteretic model, wherein the bending stiffness K is replaced in the frequency domain with the complex bending stiffness $(1+{\rm{i}}\eta )K$, where ${\rm{i}}$ is the imaginary unit [55].

Following [42], we utilize the Galerkin method to solve equation (1) with boundary conditions in equations (2a)–(2d). We use the first fundamental mode shape of a cantilever beam as a basis function, and we identify ${C}_{{\rm{A}}}$, ${C}_{{\rm{D}}}$, and M_p by fitting the experimental data in figure 5(a) using the Mathematica built-in function 'Findfit' (www.wolfram.com). Following this procedure, we find ${C}_{{\rm{A}}}=0.53$, ${C}_{{\rm{D}}}=5.91$, and ${M}_{p}=2.144\times {10}^{-5}\;{\rm{N}}\;{\rm{m}}$ with a coefficient of determination ${R}^{2}=0.9806$. To illustrate the accuracy of the model, in figure 5, we report model predictions against experimental findings.

The experimentally identified hydrodynamic coefficient is slightly lower than typical experimental observations, whereby data on ${C}_{{\rm{A}}}$ are generally found to cluster in the vicinity of one and data on ${C}_{{\rm{D}}}$ are often found in the range 10–20 [41–43, 54, 57]. Such discrepancies may be ascribed to the low aspect ratio $L/{b}_{0}$ of the IPMC, which favors 3D fluid-structure interactions [42, 43]. This hypothesis rests upon experimental observations in [42, 43], which indicate that both the hydrodynamic coefficients strongly decrease as the aspect ratio decreases, and computational results in [30], suggesting that the added mass effect is highly mitigated by 3D phenomena. Beyond the dependence on the aspect ratio, the hydrodynamic coefficients vary as functions of both the Reynolds ($\mathrm{Re}=\omega L{\delta }_{p}{\rho }_{{\rm{f}}}/\mu$), constructed by using the IPMC free length as the characteristic length in the flow, and the Keulegan–Carpenter ($\mathrm{KC}=2\pi {\delta }_{p}/{b}_{0}$) numbers, which quantify the role of viscous and convective phenomena in the flow field [41, 57–59].

The identified value of M_p is in line with results in [7]. Specifically, the value of M_p for an applied voltage of $3{\rm{V}}$ can be retrieved using the same physical and geometric parameters in [7], namely, Faraday constant ${\mathcal{F}}=96485\;{\rm{C}}\;{\mathrm{mol}}^{-1}$, gas constant ${\mathcal{R}}=8.314\;{\rm{J}}\;{\mathrm{mol}}^{-1}\;{{\rm{K}}}^{-1}$, temperature ${\mathcal{T}}=300\;{\rm{K}}$, mobile counterion concentration ${c}_{0}=1200\;\mathrm{mol}\;{{\rm{m}}}^{-3}$, volume fraction of the ionomer in the composite layers $\phi =0.5$, packing limit factor $\nu =0.25$, permittivity ratio factor ${\epsilon }^{*}=1.5$, and ionomer's semithickness $h=100\;\mu {\rm{m}}$, see [7] for detailed description. Only interfacial properties, such as ionomer's permittivity $\epsilon =2.821\times {10}^{-8}\;{\rm{F}}\;{{\rm{m}}}^{-1}$ and composite layer's thickness $d=2.85\;\mu {\rm{m}}$ are slightly modified with respect to [7] to reflect experimental observations, which are known to be highly influenced by specific fabrication conditions [48]. Specifically, the ionomer's permittivity is obtained from the measured impedance [48] and the composite layer's thickness is tuned to match the desired value of M_p.

In figures 6 and 7, we illustrate the instantaneous flow velocity and vorticity fields over one cycle of oscillation at $f=1.5\;\mathrm{Hz}$ for 8 consecutive time instants from t = 0 until $7T/8$ with step size of $T/8$, where $T=1/f$ is the oscillation period. In figure 6, the contour variable is the velocity magnitude $\parallel {\bf{u}}\parallel =\sqrt{{{\text{u}}}^{2}+{{\text{v}}}^{2}}$, normalized by the peak tip speed $\omega {\delta }_{p}$, where ${\text{u}}$ and ${\text{v}}$ are velocity component along the x and y-axes, respectively. Note that when t = 0 and $T/2$, the IPMC is not exactly parallel to the x-axis, due to a mild asymmetry of the IPMC motion caused by a slight initial waviness. In figure 7, the contour value is the normalized vorticity magnitude $\left(\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}\right)/\left({\delta }_{p}\omega /L\right)$, where positive values refer to clockwise vorticity, and negative values indicate counterclockwise vorticity.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The influence of the IPMC motion on the fluid flow is largely confined to its proximity, especially towards the tip of the IPMC, where we observe the shedding of vortices of alternating sign. In agreement with [28, 29], this results into the formation of two flow jets, as shown in figure 8(a), where we plot the mean flow field over an entire cycle computed from the full 0.66 s of recorded data (100 individual velocity fields), corresponding to 1 complete oscillation for $f=1.5\;\mathrm{Hz}$. In figure 8, we also present the mean velocity field for $f=1.1\;\mathrm{Hz}$ (figure 8(b)) and $f=2.7\;\mathrm{Hz}$ (figure 8(c)). In agreement with [28, 29], as the input frequency increases, the topology of the fluid jets changes, yet a clear dependency of the strength and angle with respect to Re is difficult to isolate. The interplay between viscous and inertial phenomena in the flow field has a central role in shaping the strength and configuration of the fluid jets originating at the IPMC tip.

Zoom In Zoom Out Reset image size

A pressure reconstruction method [36] is utilized to estimate the pressure distribution in the fluid surrounding the IPMC from the planar PIV data. Following [36], we neglect viscosity, gravity, and compressibility so that the pressure field is governed by the Poisson equation

$$\begin{eqnarray}&&\displaystyle \frac{{\partial }^{2}p(x,y,t)}{\partial {x}^{2}}+\displaystyle \frac{{\partial }^{2}p(x,y,t)}{\partial {y}^{2}}=-{\rho }_{{\rm{f}}}F(x,y,t),\end{eqnarray} \tag{ 4 }$$

where $F(x,y,t)$ is a kinematic source term written as

$$\begin{eqnarray}F(x,y,t) & = & {\left(\displaystyle \frac{\partial u(x,y,t)}{\partial x}\right)}^{2}\\ & & +2\displaystyle \frac{\partial v(x,y,t)}{\partial x}\displaystyle \frac{\partial u(x,y,t)}{\partial y}+{\left(\displaystyle \frac{\partial v(x,y,t)}{\partial y}\right)}^{2}.\end{eqnarray} \tag{ 5 }$$

Dirichlet boundary conditions are imposed along the periphery of the PIV domain, whereby we assume that the pressure vanishes away from the IPMC. In addition, neglecting the local curvature of the vibrating IPMC, a homogeneous Neumann boundary condition $\partial p/\partial x$ is imposed on the grid points close to the IPMC. The pressure field is computed using a custom Matlab script, similar to [36].

Figure 9 displays the non-dimensional pressure fields $p/\left({\rho }_{{\rm{f}}}{\delta }_{p}^{2}{\omega }^{2}\right)$ for one oscillation cycle at the same intervals and frequency with figure 6. As the IPMC vibrates, the positive pressure mostly leads the motion of the IPMC while the negative pressure follows the IPMC [53], see $3T/8$ and $7T/8$ in figure 6. As a result, the pressure difference of the fluid from the leading and trailing regions creates an opposing force to the IPMC motion [53], which is manifested through the added mass and hydrodynamic damping in equation (3). We note that the production of vorticity on the surface of the IPMC is directly proportional to the pressure gradient along the surface; thus, the concentrated regions of low and high pressure near the tip of the IPMC in figure 9 (see for example $3T/8$) correspond to the phases at which vorticity is predominantly produced and subsequently shed in figure 8 [60].

Zoom In Zoom Out Reset image size

Figure 10(a) presents the non-dimensional pressure difference across the IPMC along its length at the same intervals shown in figure 6. Specifically, the difference is calculated by subtracting the pressure on the right side from the left side of the white block in figure 9. The decrease in pressure magnitude in the vicinity of the tip of the IPMC can be ascribed to vortex shedding [30]. In agreement with [27, 30], the pressure difference is maximized in the vicinity of the IPMC tip. In figures 10(b) and (c), we display the non-dimensional pressure difference across the IPMC for $f=1.1\mathrm{Hz}$ and $f=2.7\;\mathrm{Hz}$, similar to figure 8. Numerical results in [27, 30] suggest that the specific location where the maximum is attained should be influenced by the Reynolds number. However, the rather narrow range of Reynolds numbers, from 1018 to 1844, considered in this study does not allow for resolving variations in the pressure profile as a function of Re.

Zoom In Zoom Out Reset image size

The thrust generated by the vibrating IPMC is computed using two methods: (i) integrating the reconstructed pressure on the IPMC, and (ii) performing a control volume analysis around the IPMC following [29].

Specifically, we compute the time trace of the thrust per unit width $\bar{\tau }$ by integrating the pressure obtained in section 3.3 on the outer surface of the rotating block that masks the IPMC through the Matlab built-in function 'trapz,' see figure 11. Following [27, 41], the thrust is fitted as the summation of a constant and a harmonic function with radian frequency twice the input frequency ω using the Matlab built-in function 'fit'. This behavior should be attributed to vorticity shedding at the IPMC tip, which is manifested through the generation and shedding of two coherent fluid structures per cycle [27]. Similar to numerical findings in [27], the thrust is minimized when the IPMC is in its neutral (flat) configuration and is maximized when the tip displacement is maximum [21, 23, 41].

Zoom In Zoom Out Reset image size

Alternatively, the thrust of the IPMC is computed by performing a control volume analysis in a 2D region surrounding the IPMC following [29]. In our computation, we neglect momentum flux in the out-of-plane direction and utilize only planar PIV data; the expected uncertainty of such a computation should be within 10–20% as suggested in [29]. Neglecting gravity, viscosity, and compressibility, we find

$$\begin{eqnarray}&&\tau ={\rho }_{{\rm{f}}}{\displaystyle \int }_{{\mathcal{A}}}\displaystyle \frac{\partial u}{\partial t}{\rm{d}}{\mathcal{A}}+{\rho }_{{\rm{f}}}{\displaystyle \int }_{\partial {\mathcal{A}}}u\left({{un}}_{x}+{{vn}}_{y}\right){\rm{d}}l+{\displaystyle \int }_{\partial {\mathcal{A}}}{{pn}}_{x}{\rm{d}}l,\end{eqnarray} \tag{ 6 }$$

where $\partial {\mathcal{A}}$ is the boundary of the control surface, dl is a line element of $\partial {\mathcal{A}}$, and n_x and n_y are the components of the unit outward of the control surface in the x and y directions. The control surface is approximately $78\times 88\;{\mathrm{mm}}^{2}$. The first two terms on the right hand side of equation (6) can be estimated directly from PIV data. In the evaluation of the last integral in the right hand side of equation (6), we use the pressure obtained in section 3.3 along the boundaries of the control surface. To compute the mean thrust per unit width $\bar{\tau }$ generated by the actuating IPMC, equation (6) is integrated over a period. The average of the first term on the right hand side of equation (6) is assumed to be zero by assuming that the flow field is periodic.

Following [30], we calculate the non-dimensional thrust coefficient from the mean thrust, that is,

$$\begin{eqnarray}&&{C}_{\tau }=\displaystyle \frac{\bar{\tau }}{{\rho }_{{\rm{f}}}{\omega }^{2}{\delta }_{p}^{2}L}.\end{eqnarray} \tag{ 7 }$$

Figure 12(a) presents the non-dimensional thrust coefficient with the regression model presented in [30], based on 3D direct numerical simulations. Along with predictions from [30], in figure 12(a), we also display the thrust coefficient calculated from [27] on the basis of direct 2D numerical simulations. This comparison confirms the important role of 3D phenomena and the accuracy of the regression model in [30], which can aid in the design of underwater propulsors and specifically inform the selection of the aspect ratio that is often neglected in established modeling schemes [61, 62]. We comment that significant 3D phenomena are predominantly restricted to the vicinity of IPMC tip and are more relevant for lower aspect ratios [30].

Zoom In Zoom Out Reset image size

To demonstrate the possibility of predicting IPMC thrust, we combine the distributed model in equation (1) with the regression model for the thrust coefficient illustrated in figure 12(a). Specifically, we utilize theoretical predictions in figure 5(b) to calculate the Reynolds number and use the regression model for ${C}_{\tau }$ along with equation (7) to estimate the mean thrust per unit width. Theoretical results are presented in figure 12(b), superimposed to experimental data. The model is successful in predicting the trend of variation of the mean thrust and anticipating the presence of a maximum at approximately $2.3\;\mathrm{Hz}$; the model does tend to under-predict the mean thrust for vibration frequencies in the vicinity of $2\;\mathrm{Hz}$, though the discrepancy is modest.

To compare our findings with the technical literature on thrust production of vibrating strips, in figure 13, we plot our data together with existing experimental results [21, 28, 29, 39–41]. Such data encompass IPMCs [21, 28, 29], piezoelectrics [39, 40], and servomotor-controlled passive foils [41]. The IPMC data in [21, 28, 29] share the same range of Reynolds numbers as our dataset, and consistently display similar thrust levels. The data in [21] are slightly larger than the rest; this discrepancy may be ascribed to the presence of a passive fin at the IPMC tip in [21]. Data in [39–41] show much larger thrust levels, which are in fact associated with higher Reynolds numbers. Experimental data presented therein refer to larger structures, whose characteristic length is approximately $50\;\mathrm{mm}$ [41], $90\;\mathrm{mm}$ [40], and $65\;\mathrm{mm}$ [39] as compared to $35\;\mathrm{mm}$ used in IPMCs. Servomotors and piezoelectrics afford higher tip speed compared to IPMCs, whereby data in [40, 41] encompass tip speeds as large as 0.62 and $0.18\;{\rm{m}}\;{{\rm{s}}}^{-1}$. In [39], the higher range of Reynolds number is associated with a considerably larger tip speed of $5.23\;{\rm{m}}\;{{\rm{s}}}^{-1}$, which compensates for the fact that the experiments are conducted in air.

Zoom In Zoom Out Reset image size