Collective locomotion of two-dimensional lattices of flapping plates

We study the propulsive properties of rectangular and rhombic lattices of flapping plates at O(10--100) Reynolds numbers in incompressible flow. We vary five parameters: flapping amplitude, frequency (or Reynolds number), horizontal and vertical spacings between plates, and oncoming fluid stream velocity. Lattices that are closely spaced in the streamwise direction produce intense vortex dipoles between adjacent plates. The lattices transition sharply from drag- to thrust-producing as these dipoles switch from upstream to downstream orientations at critical flow speeds. Near these transitions the flows assume a variety of periodic and nonperiodic states, with and without up-down symmetry, and multiple stable self-propelled speeds can occur. As the streamwise spacing increases, the plates may shed typical vortex wakes that impinge on downstream neighbors. With small lateral spacing, rectangular lattices yield net drag, while rhombic lattices may generate net thrust, sometimes with high efficiency. As lateral spacing increases, rectangular lattices begin to generate thrust, eventually with slightly higher efficiencies than rhombic lattices, as the two types of lattice flows converge. At Re = 70, the maximum Froude efficiencies of time-periodic lattice flows are about twice those of an isolated plate. At lower Re, the lattices' efficiency advantage increases until the isolated flapping plate no longer generates thrust. The mean input power needed to generate the lattice flows can be estimated in the limits of small and large streamwise spacings, with small-gap and Poiseuille-like flows between the plates respectively in the two cases. For both lattices, the mean input power saturates as the lateral spacing becomes large (and thrust occurs). At small lateral spacings, the rhombic lattices' input power may be much larger when the plates overlap, leading to a decrease in Froude efficiency.


I. INTRODUCTION
according to the forces on the bodies, the nonlinear dynamics are generally sensitive to initial conditions as well as the details of close interactions and/or collisions between bodies. It is very difficult to classify the whole range of possibilities in such systems. Many studies have instead focused on configurations seen in groups of biological organisms [38,45,[80][81][82][83][84][85][86][87][88][89][90][91]. Other recent studies have used machine learning to determine optimal motions of groups of swimmers [92,93]. Another large body of work concerns the use of simplified laws of interaction in place of detailed fluid dynamics, to model schools and flocks of bodies [94,95].
Following previous models [78,79,96], experiments [71,72,97], and simulations [25,[98][99][100][101][102] inspired by biology [103], we consider a particular version of the multiple flapping foil problem, with simple body geometries and kinematics, that is amenable to a wide (though by no means exhaustive) exploration of parameter space: thin plates that are oscillated vertically and moved horizontally together through a viscous fluid. The plates and motions considered here are fore-aft symmetric for simplicity; adding a pitching motion [26], an asymmetric body thickness profile [99], and/or active and/or passive deformations [73,93] can enhance the thrust generated and the propulsive efficiency.
The main quantities of interest are the time-averaged horizontal force (i.e. thrust or drag) and the input power needed to oscillate the plates vertically in the fluid. We study perhaps the most common measure of efficiency, the Froude (propeller) efficiency, a ratio of average propulsive power to average input power required to oscillate the foils [1,4].
We also study an alternative output quantity, the self-propelled speed(s) of the foils [25,27,73,93,97]. We consider a lattice of plates (rectangular or rhombic, shown in figure 1), each plate moving with the same velocity U(t) = (U, V (t)), constant in the horizontal direction, and sinusoidal in the vertical direction: with A the amplitude and f the frequency of the vertical displacement corresponding to V (t). The exponential term allows smooth start-up from rest with time constant t 0 = 0.2 (the particular value is not important, as our focus is on the eventual steady state behavior). We nondimensionalize quantities using the plate length L as the characteristic length, the flapping period 1/f as the characteristic time, and the fluid mass density ρ f as the characteristic mass density.
We solve the incompressible Navier-Stokes equations, nondimensionalized, in the rest frame of the lattice [104]: The basic dimensionless parameters are where ν is the kinematic viscosity of the fluid and L x and L y are the lattice spacings in the x and y directions, respectively. Other important dimensionless parameters, combinations of those above, are: Re is the Reynolds number based on the mean vertical velocity of the foil on each half-stroke, and is therefore a better measure of the ratio of inertial to viscous forces than Re f , which we think of as a dimensionless frequency. It is convenient for computations to nondimensionalize time by the flapping period, but the flapping frequency is one of the kinematic parameters we vary as we search for optimal flapping kinematics as well as plate spacings. Therefore, for comparison across kinematic parameters, Re U gives a more uniform measure of the horizontal speed of the foil than U L (since L and ν are considered fixed in all cases, while f varies). To find the horizontal velocities that yield efficient thrust-generating states and self-propelled states, previous work has shown that we should search in certain ranges of St (or U A , twice its reciprocal) [4,8,[29][30][31][32][33][34].
Instead of prescribing the horizontal velocity U L , one can allow it to evolve dynamically according to Newton's second law, setting the plates' rate of change of horizontal momentum equal to the horizontal component of the net fluid forces on them [25-27, 73, 100]. In this case, we have the plates' dimensionless mass M as a control parameter instead of U L . We have simulated this case, with U L (t) "free" and M fixed, and the case of fixed U L , and in both cases, periodic and nonperiodic flow dynamics can arise generically at different parameters. The coupling of body and fluid motion seems to add some additional complexity to the problem, so here we focus on the case with fixed U L , which is also the focus of most previous flapping foil studies, including those that investigated Froude efficiency [1,4,6,29]. The case with fixed U L and zero time-averaged thrust corresponds to the large-mass limit of cases with time-varying U L -those with initial conditions such that the fixed value of U L is an attracting state.
The flow starts at rest, and evolves until it converges to a periodic steady state, or remains nonperiodic up to a chosen end time of a simulation (typically t = 15 or 30). Some of these nonperiodic states may eventually converge to periodic in longer simulations. However, most have irregular oscillatory behaviors and seem likely to remain nonperiodic. These cases seem to require much longer simulations to precisely compute the long-time averages of fluid forces and input power. Thus we mostly focus on the parameters that yield a periodic state, generally those at lower Reynolds numbers, but give information about nonperiodic results in some cases.
For a plate with zero thickness in a viscous flow, the pressure and viscous shear stress diverge near the plate tips as the inverse square root of distance [105,106]. In the limit of zero plate thickness, the contribution of the pressure on the plate edges to the net horizontal force is zero. The net horizontal force on the plate is due only to the viscous shear stress on the two sides of the plate, The bracket notation denotes the jump in ∂ y u along the plate (the value at the top minus the value at the bottom).
The vertical force is due to the pressure difference across the plate: Important related quantities are the input power P in (t) and the Froude efficiency η F r : HereP in (t) is the input power nondimensionalized with ν/L 2 in place of f , for comparison across cases with different f (since L and ν are assumed fixed). The numerator and denominator of η F r both acquire factors of Re 3 f with the same change in nondimensionalization, resulting in no change for η F r .
Since ∇p in (2) is doubly periodic, it has a Fourier decomposition in which the mean (or constant part) has components we denote ∆p x /l x and ∆p y /l y , and ∇p 1 is the remainder (the mean-zero part). Thus p is decomposed into mean-zero and linear terms: p 1 is determined (up to a constant) by the incompressibility condition, ∇ · u = 0. The constant is fixed by setting p 1 to zero at an arbitrary point (e.g. the lower left corner of the flow domain). To fix the unknowns ∆p x and ∆p y , we impose a condition on the net fluid flow in the vertical and horizontal directions. We assume that the lattice of flapping plates is situated in a larger flow domain that ends at solid boundaries, where the flow is zero. We therefore assume that the spatially periodic flow approximates the flow away from the boundaries, but take the spatial average of the flow in the lattice to be zero at all times, to match that at the boundaries. The same assumption has been used in theoretical and computational studies of sedimenting suspensions at zero [107][108][109][110][111][112][113][114][115] and nonzero [116][117][118][119] Reynolds numbers, and with background turbulence [120]. In these studies the flow is typically solved in a periodic lattice or periodic cell domain, and the velocity of the sedimenting particles relative to zero-volume-flux axes is interpreted as the velocity in the physical or lab frame. In our case, the plates have zero volume, so the volume flux is that of the fluid alone. For periodic domain models of sedimentation and in the present work on flapping locomotion, there is assumed to be a transition region near the boundary where the flow deviates from spatially periodic, to obtain zero flow at the boundary. In a sedimentation simulation, Mucha et al. found that including the boundary region in the simulation had a negligible effect on particle velocity statistics far from the boundary [117].

III. NUMERICAL METHOD
We choose a flat plate geometry instead of a thin curved body (e.g. an ellipse) because it fits a periodic rectilinear grid, at the expense of creating flow singularities at the plates' edges. To study the effect of the singularity on a finitedifference solution of (2)-(3), we study a simpler model problem with the same type of singularity: potential flow past a flat plate, shown in figure 2A. The plate is the red line segment-extending along (−1/2 ≤ x ≤ 1/2 ; y = 0)-and the complex potential is w(z) = 1/4 − z 2 , with branch cut lying along the plate. We solve Laplace's equation for the stream function ψ = Im{w} in a rectangle R centered at the origin (the plate center), with lengths 3 and 2 in the x and y directions, respectively: ψ is continuous but its first derivatives diverge as inverse square roots of distance from the plate edges. Based on Stokes-flow solutions and local asymptotics of Navier-Stokes solutions [105,106], ψ has the same type of singularity as the velocity components in (2)-(3) (i.e. the viscous flows, not the potential flow defined by ψ). Both ∇ 2 ψ and the highest-order derivatives in (2)-(3), i.e. ∇ 2 u, ∇ 2 v, and ∇p, diverge as distance from the plate edges to the -3/2-power.
We will use second-order finite differences to discretize both the test problem (10) and the viscous problem (2)-(3), even though the Taylor series expansions underlying the finite difference approximations diverge at the plate edges.
Our goal with the test problem is to measure the error in a case with a simple analytical solution, given by (12) in all of R.
To mitigate the errors, we will employ nonuniform rectilinear (tensor-product) grids, and concentrate grid points near the plate edges, and along the plate surfaces, as shown in figure 2B. For the viscous problem, we use the MAC (marker-and-cell) scheme for incompressible flows [121] with the grid aligned with the plate as shown in the sample grid in figure 2C. The velocity components u and v are solved at the crosses and circles respectively, and the pressure p is solved at the triangles. The x and y locations of the symbols are either on the grid lines, or at midpoints between grid lines. For the test problem, we solve ψ on the u-grid in panel C (i.e. at the crosses).
The grid lines are defined from uniform grids using a grid-stretching parameter η. For the x-grid on the plate in panels B and C, we first define a uniform grid from -1/2 to 1/2 in X, and then the x coordinates of the grid are defined by: If the uniform grid spacing in X is ∆X, the spacing ∆x on the stretched grid x increases from approximately 1 − η at the plate edges to 1 + η at the plate center. As η increases from 0 to 1, the stretched grid transitions from uniform to highly concentrated at the plate edges. We choose a number of grid points for the plate, and then for the x-grid to the left and right of the plate in panel B, we set the number of grid points to be approximately that in the grid along the plate, scaled by the ratio of the outer region length to the plate length (unity), raised to the 1/2 power.
The functional form of the grids is the same as that in (13), with a stretching factor chosen so that the grid density is approximately continuous at the plate edges. The y grid is defined similarly to (13), given a uniform grid in Y . In the viscous computations that follow, the total numbers of grid points in x and y are similar (within a factor of 2), and in the potential flow test problem here they are equal (and denoted m). For the test problem, we solve (10) on the grid shown by crosses in figure 2B, for various values of m and η. Due to the discontinuity in flow quantities (e.g. velocity derivatives and pressure) across the plate, we use one-sided finite differences near the plate for all derivatives in both the test problem and the viscous solver, to maintain their accuracy away from the plate edges. To describe when the accuracy becomes hampered by ill-conditioning, we present the 2-norm condition number of the discrete Laplacian matrix for ψ, for various m and η in figure 2D. Each colored line plots the condition number versus m for a given η, ranging from 1 − 2 −2 (i.e. 0.75, darkest blue line) to 1 − 2 −14 (lightest yellow line), in order of increasing concentration of points near the plate edges. For each η, the condition number initially grows faster than m −2 , then asymptotes to this scaling for sufficiently large m.  We now study the effects of m and η on errors for the test problem, where the analytical solution is known. In figure 3A-C, we plot a few different measures of error in ψ. Panel A shows the infinity (sup) norm of the error in the computed ψ over the full grid relative to the exact solutionψ, given by (12) in all of R. Each colored line again corresponds to a particular η value, a few of which are labeled in the legend (the remaining values lie between these values, and are of the form 1 − 2 −k/5 for k = 10, 11, . . . , 69,70). For each line, the error initially decreases rapidly with increasing m, then much more slowly (as m −1/2 ) for further increases in m. This transition again reflects the transition in where the density of points is being increased most, near the tip at smaller m, and uniformly at large m. The singularity at the plate edges reduces the scaling from m −2 for a smooth problem with second-order finite differences to m −1/2 . However, by choosing the best η for a given m, we obtain the lower envelope of the lines in panel A, which has scaling m −3/2 , closer to the smooth case. In the viscous simulations, we need to compute the forces on the plate, which require integrating the pressure and velocity gradient, each with inverse-square-root singularities near the plate edges. For the test problem, ∂ y ψ is the analog of the velocity gradient, with the same singularity strength.
In panel B, we plot the error in its integral over the top left half of the plate, computed with second-order finite differences and then integrated using the trapezoidal rule. The exact integral is 1/2. Each curve again corresponds to a given η, which we take to 1 − 2 −16 now, closer to 1, to see the asymptotic behavior at large m better. Each curve eventually scales as m −1/2 , but by taking the minimum error over η at a given m, we can do much better. In  We take m between 256 to 512, and find that these choices are sufficient to resolve the flows and fluid forces on the plate to within a few percent in relative error for the ranges of parameters (e.g. domain size, Reynolds number, etc.) studied.
Before studying lattices of flapping plates, we examine a single flapping plate in flows of various speeds to establish a baseline of performance to which the lattice configurations can be compared. We solve the second-order finite difference discretization of (2)-(3) on the MAC grid (e.g. figure 2B and C) as a fully coupled system for u, v, and p. To simulate an isolated flapping plate in an unbounded fluid, we employ upstream, downstream, and sidewall boundary conditions in the rectangular domain, and take the boundaries sufficiently far from the body that they affect the results by less than a few percent in relative error. The upstream and downstream sides are 3 and 8 plate lengths from the plate's leading edge, and the sidewalls are 5 plate lengths from the plate. At the upstream boundary, to organize into the familiar reverse von Kármán street [8]. One panel further to the right (purple box in the third column of the top row) is St close to the self-propelled state, where F x = 0. In the last column, the flow speed is much larger and the body experiences drag although the wake still resembles a reverse von Kármán street, but with more widely spaced vortices. The second and third rows show the same sequences of flows as oncoming flow speed is increased, but at successively smaller Re f . Viscous diffusion of vorticity is more apparent, particularly in the bottom row. In the bottom two rows, the negative (blue) vortices move upwards relative to the positive (red) vortices at larger oncoming flow speeds, indicating the transition from reverse towards regular von Kármán streets [35]. For fish and birds at Re = 10 3 − 10 5 , the optimally efficient St are generally in the range 0.2-0.4 [8], while here Re = 10 − 10 2 , and the most efficient St are higher. The St that are close to optimal for Froude efficiency (green boxes) increase as Re decreases (from top to bottom rows), which is also seen in organisms as Re decreases [34].  extending from zero; this interval becomes larger as Re f increases. The curves at the lowest Re f have a U-shape to the right of zero velocity, indicating that zero velocity is an unstable equilibrium and the self-propelled state is the single stable equilibrium.
To quantify the general features of the self-propelled state, and at smaller oncoming flow speeds, the efficiencymaximizing state, we compute F x and η F r across a wide range of dimensionless frequencies (Re f ) and amplitudes (A/L). Figure 6B shows St SP S , the Strouhal numbers of the self-propelled states, where F x = 0. The numbers grow rapidly as Re f decreases to zero, and we expect divergence at a certain Re f , as there is apparently no self-propelled state below a critical Re f in experiments with rectangular plates [24] and simulations of thin ellipses [25]. For a given increases with Re f , as vortex shedding becomes more significant. The efficiency reaches a maximum of 0.06 as Re f increases to 200. Other experimental and computational studies have found the Froude efficiency is nearly unity at much higher Reynolds numbers [4,128]. Panel B also shows that at small Re f , the efficiency-maximizing A/L > 0.6, and gradually decreases to 0.2 as Re f increases to 200.
The Froude efficiency is perhaps the most common measure of efficiency in flapping-foil studies, but it is not the only way to study optimal motions. One can also consider the state that maximizes a desired output (mean thrust, showing that generally increased Re U,SP S correlates with increased P in . There is a smaller spread in the transverse direction. The Pareto frontier is the lower right boundary of the region: the set of points that maximize Re U,SP S for a given P in , or that minimize P in for a given Re U,SP S [129]. This set provides an alternative definition of maximum efficiency, that provides different optima for a range of P in , and moving in the direction transverse to the Pareto frontier, we have increasing or decreasing optimality of states. One could also replace the desired output Re U,SP S with Froude efficiency, η F r , and allow U to vary as a third input. We discuss this alternative later in the paper. other words, to change speeds, it is most efficient to vary the frequency of the plate but keep its amplitude roughly fixed. A similar trend has been observed for the tail beats of various fish species as they vary their swimming speed [130]. For an isolated flapping body, the average horizontal force is sensitive to the oncoming flow speed (i.e. minus the swimming speed) as it increases from zero. The average force is zero (or close to zero) at zero oncoming flow speed, then becomes thrust, then decreases back to zero at the self-propelled state, then becomes drag as swimming speed

V. DOUBLY PERIODIC LATTICES OF PLATES
Having established some of the main properties of an isolated flapping plate in an oncoming flow, we now discuss the much larger space of doubly periodic lattices of flapping plates. We consider two types of lattices, rectangular and rhombic (diamond), shown in figure 1, also considered by [78,80,84,86,91]. Between these two lattices is the full range of 2D (oblique) lattices. We focus only on the two endpoint lattice types (rectangular and rhombic) because the remaining parameter space is already quite large (five-dimensional). are periodic on length scales longer than a single unit cell. We focus on time-periodic dynamics for the most part, because such states are generally reached within 5-30 flapping periods. Dynamics usually appear to be nonperiodic at larger Re values, and may therefore require much longer run times to compute long-time averages with high accuracy.
We generally avoid presenting time-averaged values for these cases except where noted explicitly in the text (e.g. for the average input power). with A/L = 0.2, l y = 1, and l x ranging from 1.2 to 1.5 (labeled at right). The vertical gap is one plate length, but the horizontal gap is smaller, 0.2 to 0.5. F x initially increases with U/f A, so unlike a single flapping plate at this Re, zero velocity is a stable equilibrium here. After reaching a peak, each curve drops (sharply for l x = 1.2, then more smoothly as l x increases), and then adopts a U-shape somewhat similar to that in the isolated body case. For l x = 1.2 to 1.4, there are three zero crossings (counting U/f A = 0), corresponding to three equilibria, two stable and one unstable. Panel B shows the same data with l y increased to 1.5. The darkest blue curve (l x = 1.2) now has two sharp drops, between which the curve increases with U/f A. Near zero U/f A, the curve is dotted, indicating that the dynamics are nonperiodic in this region, so in the averages on the dotted line there is some uncertainty (that we do not quantify here). At the last of these nonperiodic cases (highest circled point), the time trace of F x (t) is shown in a small inset panel, with tick marks every flapping period along the horizontal axis. The graph of F x (t) then drops sharply to the next point, also circled, at which the time trace becomes periodic with period 1. The dynamics remain 1-periodic as F x (t) increases to the next circled point. Then F x (t) drops sharply again, to a state that is 1/2-periodic, the fourth and final circled point on this curve. The curve then increases smoothly with further increases in U/f A. In this case, the sharp drops in the curve correspond to changes in periodicity, from nonperiodic, to 1-periodic, to 1/2-periodic. We will discuss the corresponding flow structures below. The remaining curves in this panel, for l x = 1.3 to 1.5, become increasingly smooth as l x increases, eventually resembling those in panel A, but with an additional equilibrium for l x = 1.4 and 1.5; zero velocity is unstable for these cases. Panel C shows the same quantities for A/L increased to 0.5 and l y increased to 2, and a wider range of l x (labeled at right). The two inset panels show another example of the change in dynamics (from 4-periodic to 1/2-periodic) that accompany a sharp drop in F x at a particular U/f A. Panel D (l y = 3) shows three more examples of changes from 1-periodic to 1/2-periodic dynamics that occur at sharp drops in F x . Panels C and D indicate a transition with respect to l x as well. For l x near 1 (blue curves), the plates experience drag at small U/f A. At larger l x (green and yellow curves), more complex variation of F x is seen at small U/f A including thrust. By counting the numbers of zero crossings of these curves (including U/f A = 0), we obtain the number of equilibrium states, and show the totals as colored patches in the bottom row, for A/L = 0.2 (E), 0.5 (F) and 0.8 (G). In the dark blue regions U/f A = 0 is the only equilibrium, and there is no net locomotion. This is the case at smaller l y in most cases, and some larger l y values at the largest A/L (panel G). The close vertical stacking of adjacent bodies tends to suppress vortex formation and thrust generation, as we will illustrate later. In the light blue regions, there are two equilibria: U/f A = 0 is unstable and there is a stable self-propelled state, as in the isolated-body case. Examples are given by the yellow lines in panels C and D, which represent the closest approximation to the isolated body among these cases (l x and l y are largest).
However, the body might not be well approximated as isolated in some cases; there can be significant flow interactions across the periodic unit cell, particularly at A/L = 0.5 and 0.8. The yellow regions in the bottom row of panels have three or more equilibria, and these generally correspond to small l x and large l y . The interactions between adjacent bodies' edges are strongest here, and lead to a variety of flow modes (and dynamics, indicated by the insets we have discussed) that are sensitive to small changes in U/f A and the other parameters. At the largest A/L (panel G), these states are suppressed by the larger amplitude of motion, which tends to suppress interactions between vortices shed by horizontally adjacent plates.  Panels G and H also show a flow state that is up-down symmetric after a 1/2-period. Consequently, F x (t) (inset next to fourth circle in figure 10B) is 1/2-periodic-the horizontal force is the same on the up and down strokes. By contrast, the first, second, and third flow states were not up-down symmetric, and F x (t) was 1-periodic in each case.
A similar phenomenon occurs at the sudden drop in F x accompanied by the insets in figure 10C. F x (t) transitions from 4-periodic to 1/2-periodic in the insets. Flow snapshots are shown in the green box of figure 11; panels I-J for the first inset of figure 10C, and panels K-L for the second. In panels I-J, the vortex dipoles are not up-down symmetric after a 1/2-period, and the upstream member of each vortex pair is larger. In panels K-L, the dipoles are up-down symmetric, and the downstream vortices are larger. A similar phenomenon also occurs at each of the three sudden drops in F x highlighted by circles in figure 10D. The corresponding flow transitions, from up-down asymmetric to symmetric, are shown in figure 30 in appendix A. The general phenomenon then is that sudden changes from drag to thurst can occur when l x is close to 1, when the dipole jets on each half stroke switch from upstream to downstream orientations. As l x increases to larger values, the curves may becomed smoothed versions of those with sharp drops, e.g. in figure 10A. Eventually, at large enough l x , the vortex shedding pattern changes qualitatively from a dipole between adjacent leading and trailing edges, to a single dominant vortex that interacts with previously shed vortices in the wake, e.g. the reverse von Kármán street.    figure 14 where F x (t) deviates from 1-periodicity by the following measure. We compute the averages of F (t) over the last eight half-periods of V (t), during t = 11 to 15. We split the eight values into two sets of four, one for the first half-period of V (t) and the other for the second half-period. We sum the standard  deviations over the two sets and normalize by 1 4 15 11 |F x (t)|dt (an average magnitude of F x ). Where the resulting value is greater than 0.01, we plot circles in figure 14. The 0.01 threshold is somewhat arbitrary, but is chosen with certain considerations in mind. The flapping motion imposes a strong 1-periodic component in all flows, so a threshold of 0.1, say, would classify some nonperiodic F x (t) as periodic. A threshold much smaller than 0.01 would miss some F x (t) that have almost but not completely converged to periodic near t = 15. The circles occur predominantly at large l y , and more often at small l x , though they are also found at large l x . The reason is that the close spacing of plates at small l y tends to suppress complex vortical structures, and leads to more laminar, periodic flows. As we have seen, small l x leads to formation of thin dipole jets with sharp concentration of vorticity, and complicated nonperiodic dynamics can result. Many of these cases occur at small oncoming velocity (large St), reflected in the larger number of small blue circles across the panels. As we have seen for the isolated body, above a certain flow speed, a reverse von Kármán street tends to form, in many cases due to merging or other regular interactions between the leading and trailing edge vortices. The larger number of yellow circles at A/L = 0.2 (panel A) is perhaps because at a given U/f A, U/f L is smaller in panel A, so in a given flapping period, vortices do not move as far downstream relative to the plate length in this case, leading to more complicated dynamics. Also, Re is constant (20) in all three panels, so Re f decreases from left to right. To the extent that viscous regularizing effects are more controlled by Re f , they are increased moving from left to right. The corresponding data for the rhombic lattice are shown in figure 31 in the appendix, and shows similar trends but with additional nonperiodic states at smaller l y .  there they occur at large St, i.e. smaller U/f A. We have also noted a few cases of F x (t) with longer periods: 1.5 (red crosses), 2 (red circles), and 4 (red square), all of which occur at small l x . In general, the period-1 and longer-period states occur near the nonperiodic states (black squares), so the former may be intermediaries in the transition from 1/2-periodic to nonperiodic states as the parameters are varied to allow more disordered flows. of F x (t) with increasing period or non-periodic at smaller l x and larger l y is basically preserved. Compared to the rectangular lattice data, there are more 1-periodic and nonperiodic F x (t) at small l y . This is perhaps because there is more y-distance between adjacent bodies for the rhombic lattice than for a rectangular lattice at the same l y , allowing for more complex flows. One important limiting case is the large-l y limit. Here the configuration consists of well-separated 1D arrays of flapping plates. The flows around the 1D arrays are essentially the same for the rectangular and rhombic lattices. and U/f A = 3 (large enough that regular vortex arrays may be generated, and small enough that thrust may occur).
In the top row, A/L = 0.2. Moving left to right, as the gap between adjacent plates increases, the flow changes from a dipole to a vortex street. At far right, the trailing edge vortices interact more with the previously shed trailing edge vortex than with that shed at the leading edge of the adjacent plate. There is mean thrust for the flows in the first two columns of the top row, approximately zero thrust for the third, and small net drag for the fourth. In the second row, A/L = 0.5, there is again a transition from dipole jets to a vortex street, with larger vortices now. Now there is also nontrivial coupling between adjacent rows of plates. Instead of uniform flow, in the first three panels there are bands of nearly constant vorticity above and below the vortex arrays. These are shear flows with u nearly linear with respect to y and v nearly constant. They differ from those that would be seen with an isolated flapping body, and therefore, unlike the flows in the first row, they would be altered for larger l y values. In the second row, only the third column corresponds to a state of mean thrust. In the third row, A/L = 0.8, only the first column is a state of mean thrust. It is not obvious from the form of the vortex dipole why mean thrust occurs, but there are noticeable differences with the orientation of the dipole in the second panel. Again there are linear shear zones above and below the dipole jets. In the fourth column, vertical bands of same-signed vorticity form. Now the plates are coupled strongly to both vertical and horizontal neighbors through the flow.
We have shown examples of flows with l y fixed and l x from small to large. We now reverse the parameters, i.e. hold l x fixed at 2 (an intermediate value), and vary l y from small to large. Figure 19 shows examples of vorticity fields at instants of upward (and as usual, rightward) flow. Since l y may be small, there can be significant differences between rectangular and rhombic lattice flows and we show both. The differences are most pronounced at the smallest l y . In The plates now cover the full horizontal extent of the flow field, so the entire flow is forced through the small gaps between the interleaving plate edges. Consequently, much stronger vorticity is generated at the plate edges. This flow results in a net thrust force, is not temporally periodic but F x (t) has a strong 1/2-periodic component, and is also not spatially periodic on the scale of a single unit cell (i.e. containing a single plate), but of course is periodic on the the scale of a double unit cell by the definition of the periodic boundaries (the dashed box in figure 1B). This can be seen by examining the flows below the right edges of the plates. There is a small blue vortex below the top plates' lattice has a state of net thrust at both speeds (N and P). At the largest l y = 3, the rectangular (Q, S) and rhombic (R, T) flows are similar as discussed below figure 17; both pairs generate net thrust with only small differences in their magnitudes and the corresponding Froude efficiencies. To summarize, for the rectangular lattice, only net drag occurs below a moderate l y . For the rhombic lattice, by contrast, thrust can occur at very small l y , though the flows are nonperiodic and the Froude efficiency is low. We can also see that the rhombic lattice flows converge to spatially periodic on the unit cell scale as l y becomes large, which correlates with the convergence to temporal periodicity.
By computing the relative error in unit-cell-periodicity for a number of other flows (76 in all) at various parameters, we find that unit-cell-periodicity correlates with temporal periodicity in general. Both are more prevalent when Re is small, and when l x − 1 is not very small. However, there are examples of one without the other and vice versa.
Therefore, in a large lattice of plates of plates at sufficiently large Re, the flow would be expected to deviate from the periodic lattice model with unit cell periodicity. However, in most cases considered in this paper, the deviation is not very large. Below l y = 2, the solid contour lines bend sharply leftward, and efficient thrust is only obtained for l x near 1, if at all. Here the rhombic lattice is more efficient, and yields net thrust down to l y = 0.2 (yellow dashed line), as noted previously for Re = 70. Moving to panel B, there is an overall decrease in peak Froude efficiency by about a factor of 2. The solid and dashed lines deviate more at larger l y , as the now increased A/L results in more interaction between different horizontal rows of plates. For both lattice types, the peak Froude efficiency moves to much larger l x . For larger A/L, the vorticity has a larger vertical extent. For an isolated plate, figure 7B shows that U/f A should be kept roughly constant as A/L increases, for peak efficiency. This means that U/f L increases. Thus the vortices are more spread out horizontally, and it is reasonable that the plates should be more spread out to interact with vortices efficiently. The trend continues in panel C: a further reduction in peak Froude efficiency, that occurs at still larger l x .
In all cases, net thrust is obtained down to lower l y by the rhombic lattice. Another measure of performance is the maximum self-propelled speed U SP S /f A achieved by the lattice. These data are presented in figure 21 for the rectangular (top row) and rhombic (bottom row) lattices, at A/L = 0.2, 0.5, and 0.8 (left to right columns). In all panels, the highest speeds are obtained at the smallest l x , where vortex dipole formation leads to strong forces on the plates. The much smaller speeds at larger l x may be partly due to the relatively small Re, which leads to substantial diffusion of vortices as they move over a larger lattice length scale. As in figure   20, the two lattices' speeds agree better at larger l y ; the rhombic lattice achieves propulsion at smaller l y than the rectangular lattice, with moderate to large l x . An example is the local maximum in panel D at l x = 2.5 and l y = 0.5, close to the parameters for the flow in figure 19F but at lower Re. In panels B and E, there is a small band where l x = 1.3 and l y = 2.5 and 3 where the speed is greatly reduced or zero. The reason is not obvious, but may reflect a particular aspect of dipole formation at this value of Re.
When we increase Re much above 20, more flows are nonperiodic up to t = 30, and a contour plot of Froude efficiency like figure 20 would require much longer simulations to achieve reliable long-time averages. Therefore, for higher Re, we present data only for cases that meet a threshold for periodicity. We use the same criterion described in figure 14, but increase the threshold from 0.01 to 0.08. Even at this larger threshold, F x (t) is close to periodic. The larger threshold allows us to present more values, and observe the general trends more easily. These trends are not very sensitive to the particular threshold chosen (0.08). Figure 22   We have shown in figure 22 how the maximum of the Froude efficiency over U/f A varies with respect to Re, A/L, l x , and l y for flows that are periodic or almost periodic. In figure 23 we replot the same data in figure 22, adding the Re = 20 case, and showing both P in and η F r simultaneously. Thus, the four rows from top to bottom correspond to Re = 70, 40, 20, and 10, while the three columns from left to right have A/L = 0.2, 0.5, and 0.8. This figure allows us to identify which configurations are effective at different input power budgets, and describe the Pareto frontier for this set of data. In each panel we plot the data for rectangular lattices with small squares, with convex hull shown as a solid black line. The data for rhombic lattices are small diamonds with convex hull shown as a dashed black line. For each data point, the outline color gives the value of l x and the interior color gives the value of l y (listed at right). In each panel we also plot, at the panel's values of Re and A/L, P in and η F r for an isolated body. This is shown with a red cross, or if there is no thrust for the isolated body at any of the U/f A tested, the panel has a red outline. At the lowest two Re (bottom two rows), η F r for the isolated body is zero or much smaller than for the lattices. P in for the isolated body is at the lower end of the range of lattice values, near the values for the largest l x and l y . Panel J is a somewhat special case, as the large-l x lattice values did not yield thrust, but the isolated body did yield a very small net thrust-with η F r and P in both much smaller than for the lattice values shown.
For Re increased to 40 (second row from the top), there are lattice values with P in lower than that of the isolated body, and some of these-rectangular lattices with large l x and l y -yield much higher efficiency. Moving to the top row (Re = 70), the isolated body's performance is relatively improved, lying near the middle of the range of lattice values. That is, the best periodic lattice flows have efficiencies about twice that of the isolated body. However, this should be regarded as only a lower bound on the efficiency advantage of the lattice flows. Many lattice flows were not counted due to nonperiodicity. These plots also show that input power for the rhombic lattices is generally larger than for the rectangular lattices, because the rhombic lattices cover more of the flow domain horizontally, forcing the fluid through smaller constrictions between the plates. Nonetheless, the largest η F r values here, in panel A, are for rhombic lattices. We also see that P in is generally largest for smaller l x and l y (blue symbols) due to increased flow constriction, except that many of these cases are omitted because they have nonperiodic flows (especially in the top row). substituted for Froude efficiency as the measure of performance on the horizontal axes. The self-propelled speed of the isolated body is less, but sometimes not much less, than the peak speeds of the lattices in the same panel. In some cases (panels A, F, G, and J), there is not a lattice flow that has larger Re U,SP S and smaller P in simultaneously.
However, this could easily change if temporally nonperiodic lattice flows were included.

VI. MEAN INPUT POWER
The Froude efficiency and self-propelled speed depend on how the mean horizontal force varies with oncoming flow speed. For a single flapping foil, this behavior depends on the physics of vortex creation and shedding due to large amplitude flapping at a given Reynolds number. Optimal vortex creation for thrust occurs when the foil moves a certain angle of attack in the flow; this motion can be computed but is difficult to describe with a simple analytical formula [98]. The same phenomenon underlies the prevalence and optimality of St = 0.2-0.4 for flapping locomotion at high Re [31,34]. For a lattice of flapping bodies, the process is further complicated by the additional length scales of separation between bodies, and the effects of vortices colliding with downstream bodies.
One ingredient of Froude efficiency that is easier to address analytically is the time-averaged input power. We have seen in figure 9 for an isolated flapping plate that the mean input power scales as flapping amplitude and frequency Integrating the y-component of (16) over a periodic unit cell, the left side vanishes because there is zero net vertical (or horizontal) outflow. So does the viscous term, because the inverse square root singularity in the shear stress diverges too slowly to produce a vertical resultant in the limit of zero plate thickness. Hence the integral of −∂ y p over the unit cell is zero. The integral has two contributions: one from the vertical change in p across the unit cell, and the other from the jump in p across the plate, resulting in a force applied by the plate to the fluid. Hence where ∆p y is the change in pressure across the unit cell in the y direction. For the rhombic lattice, with two plates in a double unit cell (figure 1B), the integral in (18) includes both plates, and ∆p y is the change in pressure across the double unit cell in the y direction. For steady flow with dimensionless mean y-velocity 1, which relates the input power to the pressure difference across a unit cell in the y-direction. The latter can be determined analytically for the limiting cases of the flows in figure 25.
In the limit l y /(l x − 1) 1, the flow in panel A becomes Poiseuille flow in the x-interval between the plate edges (0 ≤ x ≤ 1 in panel A), and zero flow in the rest of the flow field. Poiseuille flow is a good approximation to panel A even though l y /(l x − 1) = 1/2, not very small. The lack of streamlines above and below the plates indicates slow flow there (the local density of streamlines is proportional to flow speed). In the interval 0 ≤ x ≤ 1 the flow is nearly unidirectional with a parabolic profile for v(x), and the isopressure contours are nearly equally spaced, corresponding to the constant pressure gradient of Poiseuille flow. Using Poiseuille flow to relate the pressure gradient ∆p y /l y to the net fluid flux through the unit cell, Q = 1 · l x , The last expression in (20) is written in terms of l y /(l x − 1), the parameter that sets the validity of the approximation.
The limit l y /(l x − 1) 1 is exemplified by figure 25B. This is the Stokes flow through small gaps in a periodic array of plates. The black lines are again the streamlines. They show flow converging toward the gap, and a small recirculation region centered at the midpoints of the plates. To determine ∆p y and therefore P in in (19), we can use the Stokes flow solution for a single gap in an infinite wall, derived by [105]. In the region that is much farther from the gap than its width, the pressure is approximately constant, one constant above the wall and a different constant below the wall. In figure 25B this is shown by the colored lines, isopressure contours (values at right). The solid line contours (dark blue and yellow) have pressures close to the values at distance 0.5 above and below the gaps (i.e. far from the gap). The dashed lines (dark blue and yellow) have pressures 1% above and below those of the solid lines.
Hence, in panel B above and below the "cloverleaf" regions near the gap bounded by the dashed lines, the pressure varies by less than 2%. Therefore, the pressure field is essentially the same as that far from a single gap in an infinite wall. We approximate ∆p y in (19) as the difference between the far-field pressure constants for the infinite-wall case with net flux Q through the gap, from [105]: ∆p y /Q = −32Re V /(l x − 1) 2 π . Then we have In figure 26 we plot P in for steady flows through rectangular lattices. Each row shows data for a different value of Re V , up to 6.4, close to the threshold at which the steady flow state becomes unstable for some choices of l x close to 1. The first column shows P in versus l y for different choices of l x (colored lines, values listed at right), at U/V = 0.
The second column shows the same data in rescaled variables. P in is divided by l 2 x /(l x − 1) 2 Re V , a factor that appears in both (20) and (21). On the horizontal axis, l y /(l x − 1) is used as the dependent variable. The black line shows the Poiseuille flow scaling (20), while the black cross shows the value 32/π given by (21). The agreement is almost exact.
The third and fourth columns show the same data for U/V = 2 and 8, respectively. Here the Poiseuille flow result (20) can be modified by including a dimensionless crossflow U/V in the flow equations [104], resulting in a v(x) that is linear plus exponential. In place of (20) we obtain: with Re U = U L/ν. P in in (22) tends to (20) as Re U → 0, i.e. if either Re V or U/V → 0, so the crossflow has no effect in Stokes flow (as can also be seen by linearity). So only in the three panels near the lower right corner of figure the limit l y /|l x /2 − 1| 1, using the appropriate expressions for ∆p y . The case of large l y is essentially the same as figure 25B, so for the double unit cell of the rhombic lattice, we have twice the P in of (21). We compare these results with the computed steady Navier-Stokes results in figure 27. We use only two values of U/V now, 0 and 8, and the same Re V as in figure 26.
The first column shows the unscaled data at U = 0. Unlike for the rectangular lattice, the scaling of l y changes from l y /|l x /2 − 1| to l y /(l x − 1) at small and large l y respectively, so additional columns are needed to show the data with these separate scalings. The second and fourth columns show the small-l y scalings, with black lines showing the relationships in (23). In the bottom panel of the fourth column ((Re V , U/V ) = (6.4,8)), the colored lines with l x > 2 are shifted at nonzero Re U , due to the crossflow effect discussed for the rectangular lattice (not rederived here).
There is no visible shift for the l x < 2 lines, because the effect of the crossflow cancels for the four horizontal channel flows, two in each direction, in figure 25D. The black crosses (third and fifth columns) again show the large l y /(l x − 1) values of P in . The main difference from the rectangular lattice is the large growth in P in at small l y when l x < 2.
Propulsion occurs in many of these cases, e.g. at l y = 0.5 and 0.7 (top and middle rows of figure 22). However, the peak Froude efficiency decreases noticeably when l x drops below 2 and when l y decreases in most of these cases. , above the Reynolds numbers in the present study). Nondimensionalizing (20) using ν/L in place of V , consistent with the unsteady P in results in this paper, we havẽ The small-gap flow in panel B changes to a jet flow (e.g. figure 10) as the Reynolds number rises to 20. For Re = 10-70 the flow is intermediate between viscous-dominated (resulting in (21)) and inertia-dominated. In the latter case, the momentum theorem can be used to calculateP in in the case of steady inertia-dominated (high-Re) flow.
The calculation is the same as for the steady drag on an infinite, periodically perforated plate, given in [104], section 5. 15. In the small-gap limit l y /(l x − 1) 1, the pressure drop through a unit cell of our periodic lattice becomes the same as that through the periodically perforated plate, which in our notation (and nondimensionalized by ρ f V 2 ) is The pressure is assumed to follow Bernoulli's law on the upstream sides of the plates, while on the downstream sides, where separated flow occurs, it is assumed constant and equal to that in the gap (see [104] for details). The input power follows from (19): where again,P in has been nondimensionalized using ν/L in place of V , consistent with the unsteady results in this paper. For both Stokes flow (21) and separated flow (26) with small l x − 1, P in diverges like (l x − 1) −2 , though with different prefactors.
Time-averaged input power for the rectangular lattice is plotted in figure 28 at two values of Re and A/L (labeled at left) and U/f A (labeled at top) in the ranges already considered. The unscaled P in data are shown in the first column. P in is rescaled according to the small-l y Poiseuille flow scaling (24) in the second and fourth columns, and according to the large-l y separated flow scaling (26) in the third and fifth columns, with Re (from (5)) in place of Re V . The error bars show the range of values within one standard deviation of P in , computed using the last five period-averages ofP in (t). In many cases (i.e. Re = 10, small l y ),P in (t) is periodic, so the error bar has almost zero height, and the upper and lower horizontal hash marks overlap, appearing as a single hash mark. At Re = 70 and larger l y , the vertical extent of the error bar is noticeable, and gives a measure of the nonperiodicity of the data.
Although the data are more scattered in the unsteady case, they are qualitatively similar to the steady case ( figure   26), including the shift at increasing l x in the Poiseuille flow (linear growth regime). The steady problem neglected the effects of unsteadiness (an oscillating instead of steady vertical flow) and nonlinearity in the Navier-Stokes equations.
We extended the steady model to the case of an unsteady but linearized model, for a harmonically oscillating channel flow v(x, t) = V 0 (x)e 2πit with crossflow. Analytical solutions are again possible, but more complicated than the steady case because we are back in the five-dimensional parameter space. We did not analyze the results in detail but they seemed to agree qualitatively with the small l y -results in figures 26 and 28. Despite the complications due to vortex shedding and nonperiodicity, the steady models and the fully nonlinear simulations agree qualitatively in the behavior of P in -linear growth at small l y /(l x − 1) in the second and fourth columns of figure 28, saturation at large l y /(l x − 1) in the third and fifth columns. The main discrepancy is that the Re = 10 data have a larger magnitude than the Re = 70 data. A better fit is provided by the relation P in ∼ Re 3 (typical for pressure losses at high Re, as in the third and fifth columns), rather than the ∼ Re 2 scaling for steady Poiseuille flow used in the second and fourth columns of figure 28.
The fully nonlinear P in data for the rhombic lattices are presented in figure 29. As for the rectangular lattices, the data are presented with different scalings at small l y (second and fourth columns) and large l y (third and fifth columns). In spite of the effects of nonperiodicity, there is good qualitative agreement between the unsteady and steady (figure 27) cases. In the second and fourth columns, there is a clear divergence at small l y /|l x /2 − 1| between the lines with l x > 2 (red and orange) and the remaining lines. We did not compute cases at l y /|l x /2 − 1| as small as in figure 27, however, because they are not useful for locomotion, and in the unsteady case, the parameter space is larger and each simulation requires much more computing time. Again, the Re = 10 data are generally larger than the Re = 70 data in the second and fourth columns, perhaps indicating the importance of pressure losses beyond those of Poiseuille flow at higher Re. In the third and fifth columns, the lines at various l x seem to agree reasonably well at large l y /(l x − 1).

VII. SUMMARY AND CONCLUSIONS
We have studied propulsive properties of flapping lattices of plates, at Re = 10-100, where the flows often become time-periodic within 5-30 flapping periods. This Re range is typical for submillimeter-to centimeter-scale flying and swimming organisms [96,[132][133][134][135][136], and is relevant to the increasing number of robotic flying and swimming vehicles that inhabit this size range [137][138][139][140][141]. Froude efficiency is typically much lower in this Re range than in the higher Re range typical of most fish and birds [5,6,10], so collective locomotion may be relatively more important for achieving locomotion (efficiently, or at all, at very low Re where it is no longer possible for an isolated flapping body).
In our simulations, we used a rectilinear grid with grid points concentrated near the singularities at the plates' edges. The condition number of the discrete Laplacian remains many orders of magnitude below 10 16 for the grids used in this paper, so round-off error is not a major obstacle. For a Laplace equation with similar singular behavior, we found that the method gives better than 1% accuracy in the solution and the integral of its gradient along the plate for modest mesh sizes, and converges at 3/2 order.
We first used the method to determine the propulsive properties of an isolated flapping plate in this specific context We then studied the propulsive properties of rectangular and rhombic lattices of flapping plates at low to moderate Reynolds numbers. Not surprisingly, there is a much wider range of flows than for an isolated body (which can also be regarded as the solution in the limit of large lattice spacing). When the plates are closely spaced in the streamwise direction (l x close to 1), there are sharp transitions from drag to thrust with slight increases in the oncoming flow speed (U/f A). These correspond to changes in flow modes, characterized by vortex dipoles switching from upstream to downstream directions. Some of these flows have F x (t) with periods that are various integer multiples of half the flapping period, and others are nonperiodic. In general, nonperiodicity is more common at large Re, small l x , large l y , and small U/f A, for both rectangular and rhombic lattices, with more nonperiodicity at small l y in the rhombic case.
As l x increases beyond the vicinity of 1 with large l y (akin to a 1D tandem array), the vortex dipoles transition to vortex-street wakes like those of isolated bodies, but which collide with the leading edges of plates downstream.
Varying l y from small to large with intermediate l x (fixed at 2), the flows in rectangular and rhombic lattices are initially very different, with drag-producing Poiseuille-type flows in the rectangular case, and less periodic, thrustproducing flows in the rhombic case. Interactions between vorticity at laterally-adjacent leading and trailing edges in the rhombic lattices produced thrust efficiently at small l y . As l y increased, the rectangular lattice flows eventually shed discrete vortices and produced thrust, with somewhat higher efficiencies than the rhombic lattice in many cases.
We produced maps of maximum Froude efficiency and self-propelled speeds in different portions of the four-dimensional (Re, A/L, l x , l y ) parameter space where dynamics are close to periodic in time. At fixed Re, Froude efficiency is higher at A/L = 0.2 than at 0.5 and 0.8, and the peak occurs at gradually increasing l x (and moderately large l y , ≈ 2) as A/L increases. The rhombic lattice is more efficient at small l y , and the rectangular lattice is slightly more efficient in most cases at large l y . The highest self-propelled speeds occurred at small l x and large l y for both lattice types. As Re was increased from 10 to 70, the peak Froude efficiency increased from 0.007 to 0.07. The isolated flapping body had a much lower Froude efficiency at Re = 10, eventually rising to about half that of the optimal lattices at Re = 70. However, many lattice flows are nonperiodic at Re = 70, and including these would increase the advantage over the isolated flapping body. The lattices showed similarly-sized advantages over the isolated bodies in the maximum self-propelled speeds.
The advantage in Froude efficiency occurs even with the much larger input power required for the lattices (due to the confinement of flow between the plates). In a steady model, we found that good analytical approximations to the input power could be derived for both lattice types. These relate to the flow through a small gap at small l x and moderate l y and Poiseuille flows with crossflow at small l y and moderate l x . The unsteady flows showed scaling behaviors that agreed qualitatively with those of the steady model flows.
To limit the number of parameters under consideration and the computational complexity, we have assumed all of the plates are moved together in phase (as in [71,75] but not [42], for a small group of plates). Even with this restriction, varying the spacing between the plates gives us a degree of control of the phase between shed vortices and the motions of downstream bodies with which the vortices collide. In our study, the spacing between them is not allowed to vary with time (e.g. under fluid forces [71]). Incorporating these effects would greatly expand the parameter spaces under consideration but would be natural areas for future work.

Acknowledgments
This research was supported by the NSF Mathematical Biology program under award number DMS-1811889 (S.A.).
Appendix A: Appendix Figure 30 shows examples of the vorticity fields at transitions in flows corresponding to the three sudden drops in F x in figure 10D, as U/f A is increased slightly. In each case the vortex dipoles emitted from the gaps between the plates are oriented more downstream after the transition.