3.1. Thrust generation mechanism of coordinated motion

The reciprocal motions of cylinders make the array contract and expand periodically, which accordingly induces suction and blowing flows. A single cycle of the coordinated motion is divided into contraction and expansion phases, depending on the internal area of the cylinder array. The internal area $S_{in}$ is defined as the area of an octagon in which the vertices coincide with the cylinder centres, as shown in figure 6(a). The contraction phase corresponds to a decrease in $S_{in}$ from the maximum to the minimum, and the expansion phase is the opposite. In figure 6(b), the starting times of the contraction and expansion do not coincide for various cases of phase difference. A dimensionless time $t^*$ is introduced to match the start and end of a cycle for easier comparison between the cases, as depicted in figure 6(c). The start of a cycle, $t^* = 0$, is the time when the maximum of $S_{in}$ first appears after five simulation cycles. Hereafter, only a single cycle between $t^* = 0\unicode{x2013}1$ is considered because all forces become periodic after five cycles from the beginning of the simulation; $t^*=0\unicode{x2013}0.5$ for contraction and $t^*=0.5\unicode{x2013}1.0$ for expansion.

Figure 6.

(a) Schematic of internal area $S_{in}$ of the cylinder array. Temporal change in $S_{in}$ after the fifth period with respect to (b) time $t/T$ and (c) adjusted time $t^*$: phase difference $\Delta \phi$ = (i) ${\rm \pi} /15$, (ii) ${\rm \pi} /5$ and (iii) $4{\rm \pi} /5$. Blue and red regions represent expansion and contraction phases, respectively.

The total thrust $F_T$ generated by the cylinder array is defined to be the sum of the $x$-directional forces exerted on all cylinders. Note that $F_T$ is positive along the negative $x$ axis and has a dimension of force per unit depth in a 2-D domain. The thrust coefficient $C_T$ is then calculated using the maximum speed of the cylinders $U_{max} (=2{\rm \pi} Af)$ and the array diameter $D$:

(3.1)\begin{equation} C_T = \frac{2F_T}{\rho U_{max}^2 D} = \frac{2}{\rho U_{max}^2 D} \sum_{k=1}^{N_{cyl}} ({-}F_{x,k}), \end{equation}

where $N_{cyl}$ is the number of cylinders ($N_{cyl} = 8$). To evaluate the thrust force generated per cycle, the time-averaged thrust coefficient $\bar {C}_T$ is defined by averaging $C_T$ over a cycle, $t^*=0\unicode{x2013}1$:

(3.2)\begin{equation} \bar{C}_T = \int_{0}^{1} \frac{2F_T}{\rho U_{max}^2 D} \, \mathrm{d}t^*. \end{equation}

To focus on the effect of the phase difference $\Delta \phi$, the oscillation amplitude, frequency and Reynolds number are fixed to $A=0.5$, $f=0.5$ and $Re=10$ in § 3.1 and § 3.2.

Here, $C_T$ changes dramatically over a cycle as $\Delta \phi$ varies (figure 7a). Generally, $C_T$ is negative during the contraction phase, indicating that the force exerted on the array is in the positive $x$ direction, while $C_T$ is positive in the expansion phase. This result is somewhat counterintuitive in that the total force on the cylinder array acts along the positive $x$ direction during the contraction phase, although fluid is squeezed out from the array along the positive $x$ direction by the earlier contraction of the left two cylinder sets ($i=0,1$ in figure 3a). This mechanism is in contrast to the general thrust direction produced by the contracting motion of a continuous flexible body.

Figure 7.

(a) Time history of thrust coefficient $\bar {C}_T$ of the array. (b) Time history of $-F_{x,i}$, a force acting on each cylinder in the four cylinder sets ($i = 0$–3) along the negative $x$ direction for phase difference $\Delta \phi ={\rm \pi} /5$. (c) Time-averaged thrust coefficient $\bar {C}_T$ with respect to phase difference $\Delta \phi$ ($A=0.5$, $f=0.5$ and $Re=10$).

The time history of $-F_{x,i}$, a force along the negative $x$ direction, acting on each cylinder set $i=0$–3 is presented in figure 7(b) for $\Delta \phi = {\rm \pi}/5$ which shows the largest temporal variation in $C_T$ among the three cases in figure 7(a). The sign and magnitude of $-F_{x,i}$ differ by the position of the cylinder sets ($i = 0\unicode{x2013}3$). Notably, during the contraction phase ($t^* = 0\unicode{x2013}0.5$), the magnitude of $-F_{x,i}$ for the rightmost cylinder set $i = 3$ is dominant over the other cylinder sets, causing a negative total thrust. Meanwhile, the cylinder set $i = 3$ acts positively during the expansion phase, generating the greatest thrust among the four cylinder sets. The sign of $-F_{x,i}$ for the leftmost cylinder set $i = 0$ is the opposite to that of the cylinder set $i = 3$ during most of the cycle. According to figure 7(b), the rightmost cylinder set $i = 3$ makes the largest contribution to the net thrust over a cycle among the four cylinder sets. When we conduct an additional simulation with the two cylinders of $i = 3$ removed (i.e. with only the six cylinders of $i=0\unicode{x2013}2$), $\bar {C}_T$ decreases remarkably from 0.47 to $-0.03$ under the same input parameters as those of figure 7(b), which demonstrates the critical role of the rightmost cylinder set in thrust generation.

The effect of the phase difference on the thrust generation is provided in figure 7(c). At $\Delta \phi =0$, the time-averaged $\bar {C}_T$ in (3.2) is zero due to the axisymmetric motion of the cylinder array. However, as $\Delta \phi$ increases to ${\rm \pi} /5$, $\bar {C}_T$ rises sharply. With a further increase in $\Delta \phi$, $\bar {C}_T$ decreases and finally reaches zero at $\Delta \phi ={\rm \pi}$. The phase difference that maximizes $\bar {C}_T$ is $\Delta \phi _{peak} = {\rm \pi}/5$. The phase difference in the coordinated motion is a key parameter that breaks the symmetry of the circular array and generates net thrust in one direction. Thus, understanding the effect of $\Delta \phi$ on thrust generation is essential in elucidating the physical mechanism of thrust generation based on the coordinated motion. The flow fields around the array are examined along with the cylinder motions for the optimal $\Delta \phi _{peak} (={\rm \pi} /5)$ and two other values of $\Delta \phi$ that produce smaller $\bar {C}_T$ ($\Delta \phi ={\rm \pi} /15$ and $4{\rm \pi} /5$). For these three phase differences, the corresponding time histories of $C_T$ are presented in figure 7(a).

For all three phase differences, most of the cylinders push the surrounding fluid inwards during contraction, gradually increasing pressure $p/\rho U_{max}^2$ inside the array ($t^*=0$–$0.5$, i–iv in figure 8; see also supplementary movie 1). Because the cylinder array does not form a continuous boundary, the fluid inside the array tends to escape through the gaps from the contracting array, inducing a blowing flow out of the array. The blowing flow can be identified by the outward velocity vectors in figure 8(i–iv). At the same time, the generation of the gap flow between adjacent cylinders is affected by the size of the gaps. The gap flow is relatively extensive when the gap is large, assuming that other conditions, such as the directions of cylinder motion and pressure gradient across the gap, are similar. For $\Delta \phi ={\rm \pi} /5$ and ${\rm \pi} /15$ (figure 8a,b), the gaps on the right side of the array are generally wider than those on the left side during contraction, and thus the blowing flow is dominant on the right side of the array. For the subsequent expansion phase ($t^*=0.5$–$1.0$), the expanding array sucks the surrounding fluid inwards by lowering the internal pressure (v–viii in figure 8). Opposite to the contraction phase, the gaps are wider on the left side of the array during expansion, making the suction flow prevalent on the left side of the array. As a result, for both the expansion phase and contraction phase, the overall gap flow is biased along the positive $x$ direction.

Figure 8.

Snapshots of dimensionless pressure ($p/\rho U_{max}^2$) contours and fluid velocity ($\boldsymbol {u}/U_{max}$) vectors at eight instants with equal intervals during a cycle for phase difference $\Delta \phi _{{peak}} =$ (a) ${\rm \pi} /5$, (b) ${\rm \pi} /15$ and (c) $4{\rm \pi} /5$ ($A=0.5$, $f=0.5$ and $Re=10$). The top row (i–iv) and bottom row (v–viii) show the contraction and expansion of the array, respectively. See supplementary movie 1 available at https://doi.org/10.1017/jfm.2023.1070.

Furthermore, by virtue of the high fluid viscosity and strong viscous diffusion, the boundary layer on the cylinder surface grows thick enough to overlap with the other boundary layer of the adjacent cylinder and forms a virtual fluid barrier in the gap, suppressing the penetration of fluid through the gap. This phenomenon of hydrodynamic blockage occurs between close solid structures in the low-Reynolds-number flow regime, and is a crucial mechanism for the propulsion and passive transport of tiny organisms with clustered appendages (Barta & Weihs Reference Barta and Weihs2006; Nawroth et al. Reference Nawroth, Feitl, Colin, Costello and Dabiri2010; Kolomenskiy et al. Reference Kolomenskiy, Farisenkov, Engels, Lapina, Petrov, Lehmann, Onishi, Liu and Polilov2020; Lee, Lee & Kim Reference Lee, Lee and Kim2020a; Lee et al. Reference Lee, Lee and Kim2020b; Kasoju & Santhanakrishnan Reference Kasoju and Santhanakrishnan2021; Lee & Kim Reference Lee and Kim2021). Consistent with the previous studies, it is reasonable to expect that the effect of the fluid barrier in this study is stronger for narrow gaps than for wide gaps. To confirm this, we compare the in-gap velocity profiles in figure 9, which presents the leftmost and rightmost sets of cylinders ($i=0,3$) that have the largest phase difference $3\Delta \phi$ and thus the largest difference in gap size within one cycle. In the contraction phase with outward gap flows (i–iv in figure 9), the flow through the narrow gap produced by cylinder set $i=0$ is suppressed by the strong hydrodynamic blockage and has a flat velocity profile, while the gap flow at cylinder set $i=3$ on the opposite side with a larger gap continuously develops to a larger magnitude. As the array proceeds to the expansion phase with inward gap flows (v–viii in figure 9), the narrower gap produced by cylinder set $i=3$ forms a stronger fluid barrier and has a smaller velocity magnitude than that through the gap at cylinder set $i=0$. That is, in addition to the physical size of the gaps, the hydrodynamic blockage guides the fluid towards relatively wide gaps and enhances the directionally biased distribution of the gap flow.

Figure 9.

Normalized velocity profiles $u/U_{max}$ in the positive $x$ direction for gap flows passing through the leftmost cylinder set ($i=0$, solid line) and the rightmost cylinder set ($i=3$, dashed line) during a cycle ($\Delta \phi _{{peak}} = {\rm \pi}/5$). The eight instants (i–viii) correspond to the snapshots in figure 8; panels (i–iv) for contraction and panels (v–viii) for expansion.

The velocity profile in figure 9 also shows that the fluid near the cylinder often moves in the direction opposite to that of the fluid near the centre of the gap. That is, the signs of the velocity are different between the two locations. This behaviour is because the fluid near the cylinder is affected by the thick shear layer attached to the cylinder and moves in the direction of cylinder motion, rather than following the primary direction of gap flow. For instance, when the array expands at $t^*=0.5$–0.75, the cylinder set $i=0$ and fluid near the cylinders move outwards (negative $x$ direction), while the fluid near the centre of the gap is drawn inwards (positive $x$ direction) to fill the expanding space inside the array.

The phase difference and spatial variation in gap size are responsible for the temporal trends of $C_T$ and $-F_{x,i}$ depicted in figure 7(a,b). To better understand the dynamics of cooperative array motion, the fluid forces acting on the cylinders are schematically depicted for several key frames in figure 10. As the x-directional force is mainly caused by the leftmost ($i=0$) and rightmost ($i=3$) cylinder sets, we only examine these two cylinder sets for the fluid forces. In figure 10(a) for the early stage of contraction, despite the leading inward motion of cylinder set $i=0$, the fluid forces acting on the two cylinder sets have relatively small differences in magnitude, as shown in figure 7(b). However, during the subsequent inward motion of cylinder set $i=3$ (figure 10b), cylinder set $i=3$ is strongly affected by the high internal pressure at the array centre which has been built up by the preceding inward motion of the other cylinders, and this cylinder set receives a relatively strong fluid force along the outward direction from the mutual interaction with the other cylinders. As a result, the overall magnitude of $-F_{x,i}$ acting on cylinder set $i=3$ exceeds that of cylinder set $i=0$, leading to negative $C_T$ for the entire cylinder array (figure 7a,b).

Figure 10.

Prescribed cylinder velocity vectors (grey) and fluid force vectors (red) for each cylinder during (a,b) contraction and (c–e) expansion ($\Delta \phi _{{peak}} = {\rm \pi}/5$). The fluid force vectors acting on the cylinder are only depicted for cylinder sets $i=0,3$ for simplicity.

Similar to the contraction phase, the expansion phase can be divided into a process of internal pressure change and a process in which the forces exerted on the rightmost cylinder set are affected by the internal pressure. In figure 10(c) for the early stage of expansion, the internal pressure at the array centre decreases and becomes lower than outside the array; this is clearly demonstrated in figure 8(a). The outward motion of cylinder set $i=3$ is more strongly affected by the low pressure at the array centre because cylinder set $i=3$ moves outwards after the low pressure region has formed at the array centre (figure 10d). That is, cylinder set $i=3$ is subjected to a stronger fluid force along the inward direction from the mutual interaction with the other cylinders. The overall magnitude of $-F_{x,i}$ for cylinder set $i=3$ surpasses that of cylinder set $i=0$, yielding positive $C_T$ for the cylinder array (figure 7a,b).

According to figure 8(a,b), the pressure and velocity distributions of $\Delta \phi ={\rm \pi} /15$ and ${\rm \pi} /5$ appear to be similar at each instant $t^*$, but the blowing and suction of $\Delta \phi ={\rm \pi} /15$ are not as spatially biased as those of the peak phase difference $\Delta \phi _{peak}={\rm \pi} /5$. For $\Delta \phi ={\rm \pi} /15$, the gap size difference between the left cylinder sets ($i=0,1$) and the right cylinder sets ($i=2,3$) is relatively minor throughout the cycle, and the flow is weakly biased over the cycle as a consequence. Although the internal pressure becomes strong at the late stage of contraction, the weakly biased gap flow yields smaller $\bar {C}_T$. However, when the phase difference is large (i.e. $\Delta \phi =4{\rm \pi} /5$), both the array motion and resultant gap flow are distinctly different from the previous two cases (figure 8c). Some adjacent cylinders move in opposite directions due to the large $\Delta \phi$, leading to minor changes in the internal area and pressure over a cycle. In addition, while the gap size changes gradually along the $x$ axis for the small phase differences ($\Delta \phi ={\rm \pi} /5, {\rm \pi}/15$), the gap size changes abruptly at the large phase difference ($\Delta \phi =4{\rm \pi} /5$). Due to the unorganized spatial distribution of gap sizes at the large phase difference, large gaps appear in the middle of the array and the fluid barrier is not effectively utilized. The generated flow cannot be guided into the gaps at either end, resulting in less biased flow. The minor change in the internal pressure, coupled with the unorganized gap distribution, causes a significant reduction in $\bar {C}_T$. In summary, thrust generation through coordinated motion requires a strongly biased gap flow over a cycle to produce a large time-averaged thrust, and the intensity of the biased flow is mainly determined by the non-uniformity of gap sizes, which originates from the phase difference.

3.2. Prediction of thrust generation using motion parameters

Using the relation between the formation of a spatially biased gap flow and thrust generation discussed in § 3.1, we construct a new parameter composed of several geometric variables through a simple scaling argument. This parameter captures the trend in the time-averaged thrust with respect to the phase difference. The biased gap flow can be quantified by the net volume rate of the gap flow that enters and exits the circular boundary of the array. Considering the mass conservation of incompressible fluids, the total flux of fluid expelled out or drawn in through the gaps is equal to the change in the internal area of the array. Thus, we express the total volume flux of suction and blowing as the time derivative of the internal area, $\mathrm {d} S_{in}/\mathrm {d} t$. Here, the internal area $S_{in}$ is the area of the grey octagon depicted in figure 11(a). To normalize the rate of change in the internal area, an average rate of internal area change at $\Delta \phi = 0$ is used; for simplicity, the area is assumed to be circular. For $\Delta \phi = 0$, the area changes from a maximum of ${\rm \pi} (D/2+A)^2$ to a minimum of ${\rm \pi} (D/2-A)^2$ over half a period, $T/2$. The normalized rate is then $(\mathrm {d} S_{in}/\mathrm {d} t)/(4{\rm \pi} DA/T) =\mathrm {d} S^*_{in}/\mathrm {d} t^*$, where $S_{in}^* = S/(4{\rm \pi} DA)$. For the contraction phase, in which the rate of internal area change $\mathrm {d} S^*_{in}/\mathrm {d} t^*$ is negative, an outward gap flow is generated and the total volume flux of the outward gap flow is then quantified as $-\mathrm {d} S^*_{in}/\mathrm {d} t^*$. This quantity is closely related with the temporal change in pressure inside the internal area, as exemplified in the pressure contours of figure 8.

Figure 11.

(a) Schematic showing outward unit normal vector $\boldsymbol {n}_j$ and angle $\theta _j$ at the $j$th gap. (b) Time-averaged flux coefficient $\bar \beta$ with respect to phase difference $\Delta \phi$ ($A=0.5$, $f=0.5$ and $Re=10$).

The velocity vectors in figure 8 show the amount and direction of the flow passing through the cylinder gaps, and indicate that the flux in the gap flow at each gap is affected by the size and orientation of the gap. The size of the $j$th gap depicted in figure 11(a) is denoted as $G_j(t^*)$. We assume that the ratio of the flux across the $j$th gap to the total flux across the eight gaps is approximated as the ratio of the gap size for the $j$th gap, $G_j(t^*)/G_{tot}(t^*)$, where $G_{tot}(t^*)$ is the total gap size of the eight gaps at that instant. While this assumption is acceptable at the low Reynolds number with strong hydrodynamic blockage, it is invalid at the high Reynolds number, which is discussed in § 3.4. The dimensionless flux of each gap flow is then expressed as a product of the rate of change in the internal area and the gap ratio: $-\mathrm {d} S^*_{in}(t^*)/\mathrm {d} t^* \times G_j(t^*)/G_{tot}(t^*)$. To establish a relation between the biased gap flow and thrust generation, the degree of bias in the generated flow is quantified by the $x$-directional component of the flux in the gap flow. Considering the outward-facing normal vector $\boldsymbol {n}_j(t^*)$ of the $j$th gap, the $x$ component of the gap is $G_j(t^*) \cos \theta _j(t^*)$, where $\theta _j$ is the angle between the vector $\boldsymbol {n}_j$ and the $x$ axis. The $x$-directional flux of the flow is obtained by replacing $G_j(t^*)$ with $G_j(t^*) \cos \theta _j(t^*)$. A flux coefficient $\beta (t^*)$ is then given by summing the $x$-directional fluxes through all gaps:

(3.3)\begin{equation} \beta(t^*)=\sum_{j = 1}^{N_{gap}}-\frac{\mathrm{d}S^*_{in}(t^*)}{\mathrm{d}t^*}\frac{G_j(t^*) \cos \theta_j(t^*)}{G_{tot}(t^*)}, \end{equation}

where $N_{gap}$ is the total number of gaps ($N_{gap} = 8$). The sign and magnitude of $\beta (t^*)$ denote the direction and strength of the biased gap flow: positive when the net gap flow is along the positive $x$ direction and negative when the net gap flow is along the negative $x$ direction. Accordingly, the time-averaged value of $\beta$ over a cycle (i.e. time-averaged flux along the $x$ direction) is expected to be correlated with the time-averaged thrust coefficient $\bar {C}_T$ in (3.2).

To evaluate the time-averaged extent of the biased flow over a cycle, the time-averaged flux coefficient $\bar \beta$ is considered:

(3.4)\begin{equation} \bar\beta = \int_{0}^{1}\beta(t^*)\,\mathrm{d}t^*. \end{equation}

Indeed, $\bar \beta$ in figure 11(b) exhibits a very similar tendency to the time-averaged thrust coefficient $\bar {C}_T$ in figure 7(c). Near phase differences of $\Delta \phi =0$ and ${\rm \pi}$, both $\bar {C}_T$ and $\bar \beta$ are close to zero, and more importantly, the value of $\Delta \phi$ that maximizes $\bar \beta$ is ${\rm \pi} /5$, matching the optimal phase difference for $\bar {C}_T$. Negative values of $\bar \beta$ are also observed near $\Delta \phi =3{\rm \pi} /5$. According to the definition of $\bar \beta$, its negative value means that the array has wider gaps facing the negative $x$ direction during contraction, while during the expansion, wider gaps face the positive $x$ direction. This array motion is not ideal for thrust generation along the negative $x$ axis and, indeed, the magnitude of $\bar {C}_T$ corresponding to the negative $\bar \beta$ is very small albeit positive. Interestingly, with only time-dependent geometric variables ($S_{in}$, $G_{tot}$, $G_j$ and $\theta _j$) in (3.3) and (3.4), which can be calculated from the input parameters without numerical simulations, the trend in the time-averaged thrust generation for the complex coordinated motion of the cylinder array can be successfully predicted.

On the one hand, the time-dependent $\beta$ in (3.3) is insufficient to predict the instantaneous thrust generation of the cylinder array. For instance, $\beta$ is positive for both contraction and expansion phases, while the signs of $C_T$ in each phase are opposite to each other; see figure 17 in the Appendix. For the instantaneous thrust, the temporal change of the flow inside the array (inside the internal area depicted in figure 6) should be considered as well as the flow through the gaps (through the boundary of the internal area). On the other hand, $\beta$ does not account for the time-varying flow within the array, making it unable to predict the instantaneous thrust. However, because of the periodicity of cylinder motion and generated flow, the change of the flow within the array over a period is negligible. For this reason, while $\beta$ is not in accordance with $C_T$, its cycle-averaged value $\bar {\beta }$ is found to be correlated with the cycle-averaged thrust $\bar {C}_T$.

3.3. Effects of oscillation amplitude and frequency

In addition to the phase difference effects addressed in § 3.1, we analyse how thrust generation is affected by the oscillation amplitude $A$ and frequency $f$ while maintaining the fluid kinematic viscosity unchanged. The amplitude ranges between $A = 0.1$ and 0.5 to prevent any collisions between the cylinders, which may occur for amplitudes exceeding $A = 0.5$. The frequency varies from $f = 0.25$ to 1.0.

When the amplitude varies in the range of $A = 0.1$–0.5 for a fixed frequency $f = 0.5$, the time-averaged thrust coefficient $\bar {C}_T$ maintains a consistent trend with respect to the phase difference, having the maximal value near the phase difference $\Delta \phi = {\rm \pi}/5$, but it gradually increases with $A$ near $\Delta \phi = {\rm \pi}/5$; refer to figure 12(a-i). Specifically, for $\Delta \phi = {\rm \pi}/5$, $\bar {C}_T$ increases from 0.29 to 0.47 as $A$ increases from 0.1 to 0.5. The change in $\bar {C}_T$ is relatively minor, compared with the five times increase in $A$. This thrust increase is not solely due to the greater velocity of cylinders, which results in stronger flow velocity and pressure (in dimensional forms). The change in spacing of an array also contributes to thrust generation. In the radial oscillations of cylinders, the differences in gap sizes at a given time magnify as $A$ increases, which is crucial for enhancing the degree of bias in the generated gap flow. The effects of $A$ on the bias of the gap flow are described in figure 13(a–c) and supplementary movie 2. For a large amplitude ($A = 0.5$), the differences in the gap sizes become more pronounced in the middle of contraction or expansion. Thus, the gap flow is stronger through the wider gap between the cylinder set $i = 3$ at $t^* = 0.25$ and the wider gap between the cylinder set $i = 0$ at $t^* = 0.75$, yielding notable spatial variations in the gap flows. Note that the dimensionless pressure $p/\rho U_{max}^2$ inside the array is the most significant for $A = 0.1$ (figure 13a–c). However, strong dimensionless internal pressure itself does not guarantee the greater thrust coefficient, and the bias of the gap flow should be also considered.

Figure 12.

(a) Time-averaged thrust coefficient $\bar {C}_T$ and (b) time-averaged flux coefficient $\bar {\beta }$ versus phase difference $\Delta \phi$ for (i) five oscillation amplitudes $A$ ($\,f = 0.5$) and (ii) five oscillation frequencies $f$ ($A = 0.5$).

Figure 13.

Snapshots of dimensionless pressure ($p/\rho U_{max}^2$) contours and fluid velocity vectors ($\boldsymbol {u}/U_{max}$) in (i) the middle of contraction ($t^* = 0.25$) and (ii) the middle of expansion ($t^* = 0.75$) for (a–c) three amplitudes $A=0.1$, 0.3 and 0.5 ($\,f = 0.5$) and (d–f) three frequencies $f=0.25$, 0.5 and 1.0 ($A = 0.5$). The phase difference is $\Delta \phi ={\rm \pi} /5$. See supplementary movie 2 for panels (a–c) and movie 3 for panels (d–f).

The time-averaged flux coefficient $\bar \beta$ behaves similarly to $\bar {C}_T$ in that it increases monotonically with $A$ near $\Delta \phi = {\rm \pi}/5$. In contrast to $\bar {C}_T$, $\bar \beta$ is simply proportional to $A$ at any phase difference because the rate of change in the internal area and the gap size ratios, which constitute the flow coefficient in (3.3), are directly affected by the change in $A$. Although $\bar \beta$ does not accurately quantify the extent of thrust change by the varying amplitude, it suggests that a cylinder array with a larger amplitude produces a more biased flow.

As the frequency $f$ increases for a fixed amplitude $A = 0.5$, $\bar {C}_T$ rather decreases despite the increasing cylinder speed, according to figure 12(a-ii). With the fourfold increase in $f$ from 0.25 to 1.0, $\bar {C}_T$ decreases from 0.63 to 0.41. For the large $f$, the faster movement of the cylinder reduces the thickness of the boundary layer on the cylinder surface due to the weaker viscous diffusion, which reduces the strength of the virtual fluid barrier in narrow gaps. Thus, the flow generated by the array is not well guided to wide gaps for the higher $f$, and the gap flow through the narrowest gap becomes stronger than with the lower $f$. In figure 13(d–f) and supplementary movie 3, this result is clearly depicted by the sizes of arrows representing the dimensionless fluid velocity $\boldsymbol {u}/U_{max}$ in the narrow gaps. The weakened hydrodynamic blockage lowers the degree of bias in the gap flow and consequently decreases $\bar {C}_T$. It is noteworthy that $\bar {C}_T$ increases notably as $f$ changes from 0.375 to 0.25, while the variation in $\bar {C}_T$ is relatively minor in the range $f = 0.375$–1.0. The viscous diffusion length $\delta$ (i.e. boundary layer thickness) from a solid object generally scales as $\delta \sim (\nu /f)^{1/2}$ for oscillating flow and this scaling indicates non-uniform variation in $\delta$ with respect to $f$. For the low frequency of $f = 0.25$, the significant growth of viscous diffusion length eventually leads to the notable increase in $\bar {C}_T$.

As $f$ varies for a given $A$, the configuration of the cylindrical array remains identical at a given dimensionless time $t^*$. The rate of change of the internal area in the dimensionless form, $\mathrm {d}S^*_{in}/\mathrm {d}t^*$, as well as the gap ratio $G_j(t^*)/G_{tot}(t^*)$ of each gap do not change at the given $t^*$. Consequently, the flux coefficient $\beta$ in (3.3) and the time-averaged flux coefficient $\bar {\beta }$ in (3.4) are independent of $f$; see figure 12(b-ii). This is because the flux coefficient is determined by the geometric characteristics of the array and does not account for the viscous effects of the fluid, which is responsible for the variation in $\bar {C}_T$ by $f$ (figure 12a-ii).

3.4. Thrust generation under different flow regimes

As the virtual fluid barrier that contributes to generating thrust is known to be effective in the low-$Re$ regime, we investigate the effect of viscosity on thrust generation. In addition to $Re=10$ considered thus far, $Re=5$ and $100$ are examined as the representatives of higher and lower viscosities, respectively. The motion parameters remain the same as in § 3.1 ($A=0.5$ and $f=0.5$). No adjustments to the numerical settings were necessary for achieving convergence in the simulations with the changed Reynolds numbers; only a few additional iterations were required within each time step.

As the Reynolds number decreases from 10 to 5, the trend of thrust generation with respect to the phase difference remains similar, having the optimal phase difference at $\Delta \phi = {\rm \pi}/5$, as shown in figure 14. The enhancement of thrust is also observed across most phase differences. This result is expected from the thrust-generating mechanism explained in § 3.1. Although the virtual fluid barrier is used to create a biased flow at both Reynolds numbers, flow penetrating through narrow gaps becomes more suppressed as viscosity increases. Instead, a significant amount of fluid prohibited by the narrow gaps passes through the wide gaps, thereby increasing the degree of bias in the gap flow and producing stronger thrust. Since the same hydrodynamic blockage mechanism applies for both $Re = 10$ and 5, hereafter, we will focus on examining how the disappearance of hydrodynamic blockage affects flow distribution and thrust by comparing $Re = 10$ and 100.

Figure 14.

Time-averaged thrust coefficient $\bar {C}_T$ for $Re=5, 10$ and $100$ ($A=0.5$ and $f=0.5$).

When $Re$ increases from 10 to 100, $\bar {C}_T$ decreases at all phase differences. Obviously, the high Reynolds number ($Re=100$) produces a markedly different trend in $\bar {C}_T$ with respect to the phase difference $\Delta \phi$ compared with the low Reynolds numbers ($Re=5$ and 10) (figure 14). The value of $\bar {C}_T$ reaches its maximum near $\Delta \phi ={\rm \pi} /15$ and drops sharply until $\Delta \phi ={\rm \pi} /5$, remaining low up to $\Delta \phi ={\rm \pi}$. With $\Delta \phi ={\rm \pi} /5$, in particular, $\bar {C}_T$ decreases from 0.47 to 0.046 (figure 14). For a single isolated cylinder model, as $Re$ increases from 10 to 100, the shear stress on the cylinder wall and the resultant viscous force acting on the cylinder are reduced because of the tenfold decrease in fluid viscosity. However, for our cylinder array, the reduction in the magnitude of the viscous force due to the lower fluid viscosity is not a primary reason for the remarkable decline in $\bar {C}_T$. Similar to increasing the frequency (§ 3.3, figure 12a-ii), the weaker thrust generation for $Re=100$ can be attributed to the decline in the hydrodynamic blockage; the increase in the frequency has the same effect as the increase in the Reynolds number. The viscous diffusion length of the boundary layer scales with the square root of the fluid kinematic viscosity (Schlichting & Gersten Reference Schlichting and Gersten2000; Nawroth et al. Reference Nawroth, Feitl, Colin, Costello and Dabiri2010). With a high Reynolds number ($Re = 100$) and one-tenth of the kinematic viscosity, the boundary layer becomes notably thinner due to reduced viscous diffusion, causing the virtual fluid barrier in the gaps to almost disappear. The significantly weakened viscous diffusion and fluid barrier effect make the temporal variation in the internal pressure less distinct, and the gap flows generated through the narrow gaps are as strong as those through the wider gaps. This results in a decrease in the degree of bias in the gap flow; compare figure 15(a) with figure 15(b) and see supplementary movie 4.

Figure 15.

Snapshots of dimensionless pressure contours ($p/\rho U_{max}^2$) and fluid velocity vectors ($\boldsymbol {u}/U_{max}$) in (i) the middle of contraction ($t^* = 0.25$) and (ii) the middle of expansion ($t^* = 0.75$) for (a) $Re = 10$ and (b) $Re = 100$ ($A=0.5$, $f=0.5$ and $\Delta \phi ={\rm \pi} /5$). See supplementary movie 4.

The distinct trend of $\bar {C}_T$ for $Re=100$ in figure 14 cannot be accurately captured using $\bar {\beta }$. This limitation arises because $\beta$ is based on the assumption that the virtual fluid barrier is sufficiently effective and does not directly reflect the variations in kinematic viscosity. The decrease in viscosity and weakened hydrodynamic blockage cause substantial changes in gap flow, and the basic assumption of $\beta$ that the flux through the gap is correlated with the gap size is no longer valid.

The weak fluid barrier for $Re=100$ is identifiable from the large velocity vectors in the narrow gaps in figure 15(b), and is more apparent in the $x$-directional velocity profiles of gap flows presented in figure 16 for $\Delta \phi ={\rm \pi} /5$. As $Re$ increases from 10 to 100, the gap flow through the narrower gaps is enhanced in both the contraction and expansion phases. In the middle of contraction, at $t^* = 0.25$ (figure 16i), the gap at cylinder set $i = 0$ is narrower than that at $i =3$. By increasing $Re$ to 100, the magnitude of $u/U_{max}$ increases by a factor of 3.3 at the midpoint of the narrower gap. Similarly, in the middle of expansion, at $t^* = 0.75$ (figure 16ii), the gap at cylinder set $i=3$ is narrower than that of $i=0$, and the magnitude of $u/U_{max}$ increases by a factor of 2.5 at the midpoint of the narrower gap. In return, for both the contraction and expansion phases, the lower volume of flow is guided to the wider gap in the opposite side, causing the $u$-velocity profile to have a notable deficit near the midpoint of the wider gap. The magnitude of $u/U_{max}$ at the midpoint of the wider gap decreases from 2.84 to 0.68 at $t^* = 0.25$ and from 2.34 to 1.13 at $t^*=0.75$. As a result of the weaker hydrodynamic blockage effect, the magnitude of the instantaneous force acting on each cylinder decreases, although the shape of the temporal force profile over a cycle is retained. In summary, the lower level of viscous diffusion of the high-Reynolds-number flow hinders the fluid inside the array from exiting through the wider gap, eventually annihilating the bias of the gap flow and decreasing the time-averaged thrust.

Figure 16.

Normalized velocity profiles $u/U_{max}$ in the positive $x$ direction for gap flows passing through the leftmost cylinder set ($i=0$, solid line) and the rightmost cylinder set ($i=3$, dashed line) during a cycle; $Re = 10$ and 100. The velocity profiles are in (i) the middle of contraction ($t^* = 0.25$) and (ii) the middle of expansion ($t^* = 0.75$).