Mixing in arrays of villi-like actuators

Motile villi—small fingerlike projections in mammalian small intestines, as shown in Fig. 1—are known to exhibit a variety of contractile movements driven by smooth muscle fibers originating from the mucosa under neural-enteric control,¹ coordinated to enhance digestion through active motions. Designs, which can replicate the essential flow features generated by villi, have the potential to support the design of a range of bioinspired mixers and sensors, in a similar manner to arrays of artificial cilia^2–4 but with a relatively simpler actuator to design and control.

For example, microbial fuel cells (MFCs), a type of bio-electro-chemical reactor,^5,6 generate electrical energy from organic material within electrode-containing chambers. Although MFCs are a promising technology for clean waste management and energy production, the issue of sedimentation of particulates prevents continuous operation from being realized.^7,8 To alleviate these issues, artificial villi or simplified actuators could be embedded onto the surface of MFCs anodes for particulate sensing, dispersion, and increased power output gains via an increased surface area.

As a second example, artificial villi could also be used for propulsion in the design of novel bioinspired swimmers, by mounting arrays of motile villus-shaped protrusions onto a robot. Since artificial applications have flexibility of choice of inertial scaling through choice of forcing, knowledge of effective control strategies provides design cues for the development of bioinspired artificial prototypes at a range of length scales in a wide variety of environments. Often, the use of villi-like protrusions or cilia is associated with micro-swimmers due to the length scales of most biological cilia.⁹ Due to this fact, many studies modeling cilia,^10–12 both experimentally and numerically, are in the viscous regime, where $Re→0$⁠. However, natural systems (such as the propulsion of ctenophores) can employ cilia with lengths of the order of millimeters, with Reynolds numbers approaching 10³ (see, e.g., Refs. 13 and 14), which require modeling that incorporates inertial effects in the flow.^15,16 At even higher Reynolds numbers where inertial scales and associated momentum advection become dominant [between $Re=O(103)$ and $Re=O(105)$], two-dimensional^17,18 and experimental¹⁹ studies have shown that schools of fish that symmetrically beat tails in anti-phase gain significant energy savings compared with in-phase. Considering these results together, there is an implication that the benefits of phased oscillating motions apply to flows over a wide range of length scales and therefore across both viscous- and inertia-dominated situations. The behavior of these flows can be demonstrated in simplified two-dimensional computational models.

Artificial villi also have the potential to improve the performance of chemical sensors through controlled circulation. This is because chemical sensors often rely on close proximity of the target analyte in order to detect a measurable physical change (e.g., conductivity).²⁰ By embedding chemical sensors on the tip of artificial villi and varying the manner in which fluid is circulated, there is the possibility to control and enlarge both the sensitivity and measurement domain. Further, by embedding multiple sensors along the villus shaft, it would be possible to use differential sensing to map environmental conditions in greater detail.

While the topology of static obstructions on the walls of fluidic environments has been shown to affect the resulting flow and transport properties,²¹ knowledge of how coordinated movements affect fluidic conditions is required to develop bioinspired artificial applications effectively. Recent in silico²² and ex vivo²³ studies have investigated the hydrodynamic response of groups of villi. For example, in Refs. 24 and 25, groups of oscillating villi coupled to lid-driven cavity flow are investigated using a multi-scale Lattice-Boltzmann model in two and three dimensions. Although there is no specific biological justification for the groupings, the work illustrates how villus motility may influence transport at length scales as long as the scale of the gut. In a separate simulation study,²⁶ it has been shown that large-scale gut wall contractions lead to “microfolds” that actively oscillate the intervillus spacing, leading to the formation of vortices. Although these studies propose different muscle groups for the control of peripheral mixing, both show that coordinated changes in intervillus volume can significantly augment flow conditions. Thus, strategic control of motile villus-shaped protrusions can lead to the formation of effective peripheral mixing layers.

Significant mixing and energy-saving benefits can also be associated with phased oscillating motions across hydrodynamics scales. In the non-inertial regime [$Re=O(10−3$⁠)], numerous authors have investigated mixing by beating ciliated structures across a surface in two^27,28 and three dimensions.^29–32 These studies show that asymmetrical beating patterns are required to induce chaotic mixing and out-of-phase oscillations significantly shorten mixing timescales.

To extend more generally to artificial applications, we seek to characterize the hydrodynamics of villi-like structures throughout a range of flow regimes and provide a mapping for behaviors to be expected near a small intestine-like wall. At low inertial scales (i.e., ignoring gut-wall contractions), dominating viscous effects require symmetry-breaking motions³³ to induce mixing on any reasonable timescale. Conversely, although high inertial scales lead to turbulent effects, which necessarily induce mixing,^34–36 flows in which inertial effects are significant but not so high that hydrodynamic turbulence is present might be sensitive to the specific motility patterns employed.

In this paper, we perform a computational fluid dynamics numerical simulation of the two-dimensional hydrodynamic response to the motion of linear arrays of villi-like structures. We use a Reynolds number range that starts at values that are finite but achievable in some biological systems^13,14 (⁠$Re=10$⁠) and ends at values that are far beyond those that are feasible biologically, but that can be achieved with artificial actuators (⁠$Re=6000$⁠). We explore how coordinated in-phase and anti-phase motions affect mixing within a peripheral region near the villi, across this large range of Re. We restrict the phase between the motions of neighboring villi to the purely in-phase and anti-phase cases for two reasons: first, these are the two limiting cases of minimum and maximum change in intervillus volume over the an oscillation cycle; second, both of these cases are spatio-temporally symmetric with a zero mean, and so neither induces any first-order transport; that is, the motion of fluid over a time longer than the oscillation period is due to nonlinear effects induced by inertia. These two facts allow us to quantify the transport and mixing induced by second-order effects.

Our results indicate that larger changes in intervillus volume greatly enhance the underlying mixing mechanics generating the peripheral mixing layer. We believe our work supports the development of artificial mixers and sensors inspired by active villi. The restriction of our simulations to two dimensions means that some three-dimensional effects that may significantly impact transport—such as streamwise tip vortex production from pillar-like protrusions—are not captured. The results presented here may not be directly relevant in a flow where there is large variation in the third or spanwise dimension (either geometrically via the shape of the villi, or kinematically via the relative motion of the villi). However, the framework presented here, where mixing is enhanced by increasing the intervillus volume change, may be a useful one to explore in a three-dimensional scenario. Further, the results may be directly applicable to flows with more paddle-like actuators, or where sheets of pillar-like villi have motions that are synchronized in the spanwise direction.

This paper is organized as follows. We begin, in Sec. II, with a summary of the mathematical model, a description of the numerical method, and a review of relevant supporting analysis concepts; in particular, the mixing mechanisms typically observed in transitional flows. We then present the results in Sec. III, beginning in Sec. III A by exploring the effects of active intervillus volume change, via simulations of a single pair of villi oscillating both in-phase and anti-phase. Following that, in Sec. III B, we scale the problem to arrays of 15 villi revealing new coupling effects that enhance mixing. Note that our aim is not to optimize outputs as array size varies, but rather to identify the qualitative flow structures that are (or are not) common to the two different scenarios. We discuss the results in Sec. IV and end with concluding remarks in Sec. V.

In this section, we describe the mathematical model used to simulate fluid flow consisting of a linear array of oscillating villi-like actuators in a two-dimensional viscous fluid. To aid discussion in the following sections, we follow with a summary of the basic flow structures and distinct flow regimes that enable mixing and transport, which are observed in our results.

Two-dimensional simulations of the fluid flow generated by villi-like structures, rotating about a point near their base in an oscillatory manner, were conducted using a sharp-interface immersed boundary method. A detailed description of the code, and its validation for use in coupled fluid-structure interaction problems such as that studied here, is provided by.³⁷ The pertinent details are recalled here.

The code solves the incompressible Navier–Stokes equations

where u is the velocity field, t is time, ρ is the fluid density, p is the pressure field, ν is the fluid kinematic viscosity, and $Ab$ models the presence of the immersed boundaries, essentially coupling the fluid and solid phases and enforcing the boundary conditions for the fluid at the villi-like structure surface.

The equations are spatially discretized on a simple Cartesian mesh using a second-order central-difference scheme. The use of a Cartesian, non-body-conforming mesh simplifies the mesh construction and also allows for simple domain decomposition for distributed-memory parallel computation.

The velocity and pressure field values are solved for at the nodes this co-located mesh. Note that the “face” values of the velocity field (the values at the mid-point between the nodes) was also calculated and used to form the right-hand side of the Poisson equation for the pressure (described below in Sec. II A 2) to avoid odd-even decoupling or “checkerboarding” of the pressure.

Temporal integration of the resulting difference equations is performed using a two-way time-splitting scheme. This splits the time integration into two sub-equations, the first integrating the velocity at the beginning of the time step to an intermediate time by only considering the advection and diffusion terms {the first and third terms on the right-hand side of the first equation in equation [Eq. (1)]}. The time-splitting scheme follows the process outlined in Ref. 38.

This intermediate velocity field is then updated by imposing the motion of the villi-like structures (the immersed boundaries) over the time step {the last term on the right-hand side of the first equation of equation [Eq. (1)]}. This requires a process of updating values that cross the immersed boundary—those going from the fluid into the solid (termed “ghost” points), and vice versa (termed “fresh” points). Ghost and fresh points are assigned a velocity value via interpolation using a stencil that incorporated the immediately neighboring points and the value of the boundary condition at the closest surface point, a process that needs to be solved iteratively as many of the fresh and ghost points neighbor each other and their values are therefore coupled. The immersed boundary is represented by linear finite elements that maintain a second-order accuracy. A detailed description of this identification and interpolation process is provided in Ref. 37.

We note here that this direct use of the value of the boundary condition in the interpolation for ghost and fresh point values results in a sharp representation of the boundary. Forces are not distributed over neighboring points to weakly enforce the boundary conditions; the boundary condition is imposed exactly up to the order of the discretization (which here is second order). A sharp-interface representation is important in flows such as those studied here where the production of vorticity and the roll-up of the boundary layer into vortices dictates the dynamics.

A second sub-step equation is then formed to integrate the velocity field from this intermediate time to the end of the time step, by integrating the pressure term. First, the divergence of this equation is taken, and the divergence-free condition from continuity [the second equation in equation Eq. (1)] is imposed at the end of the step, resulting in a Poisson equation for the pressure correction (the change in pressure over the time step) with a right-hand side in terms of the intermediate velocity field. This right-hand side is formed using the face values of this field (the values interpolated from the nodes to the mid-points between nodes) to ensure odd and even points are coupled.

A further correction to the right-hand side is imposed to improve the mass conservation. Those cells of the underlying mesh that are “cut” by the boundary are identified and classified as small (where less than half the cell volume remains in the fluid) and large (where half or more of the cell volume remains in the fluid). For large cells, there is no treatment. However, for small cells, their right-hand side contribution is distributed to surrounding cells in the domain. This ensures that discrete volume is not lost as cells cross from the fluid into the solid, and significantly improves mass conservation. The scheme, referred to as a “cut-cell” technique, has previously been described, tested, and validated.³⁹ The resulting Poisson equation is then solved for the pressure correction using a geometric multigrid method, typically to a tolerance such that residuals are $<10−6$ of the maximum pressure in the domain.

Finally, this pressure correction is added to the pressure, and the velocity field is integrated from the intermediate time to the end of the time step.

The code has previously being employed for rotating⁴⁰ and multiple-body fluid–structure interaction problems,^41,42 and stability and structural sensitivity studies.⁴³ 

At the structure surface, a Dirichlet condition is imposed for the velocity from the prescribed motion of the structure. For the pressure, a zero-normal-gradient condition is imposed. More complex conditions for the pressure gradient that are derived from the Navier–Stokes equations⁴⁴ have been tested but provide effectively no change to the flow, and the current condition has been validated previously against other codes in similar fluid-structure flows.³⁷ These boundary conditions are incorporated into the immersed boundary interpolation described above in Sec. II A 2.

At the domain boundaries, a zero-normal-gradient condition is imposed for the velocity component tangential to the boundary, while the normal component is set to zero (i.e., a zero-penetration condition). This zero-stress, or Robin, boundary condition is also applied under the villi-like structures, effectively removing the boundary layer by mimicking a free surface. The domain boundary conditions are imposed at face locations, in between two mesh nodes. This leaves a layer of ghost nodes that are outside the domain. Boundary conditions are imposed by extrapolation (of the same order as the underlying scheme) for values on the ghost nodes from the neighboring points in the fluid domain and the boundary condition.

For the simulations of this study, a rectangular fluid domain of length 30 L and height 10 L—where L is the length of the villi (see below)—was discretized using a 2048 × 1024 grid, with variable mesh spacing used to concentrate points in the vicinity of the villi where high flow gradients are expected, as shown in Fig. 2. No background flow was imposed, with the flow being quiescent other than the motion induced by the motion of the villi. We reiterate that the motion of the villi is prescribed and is not a function of the fluid motion—loosely mimicking biological villi whose motion may be controlled by muscular actuation. Villi were represented as rigid rectangles of length L and width $0.25L$⁠, with semi-circular end caps, which then rotated about the center of curvature of the bottom cap. Villi were placed such that a gap of $0.1L$ was maintained between the villi and the bottom boundary of the flow domain. While this allowed a small amount of fluid to pass between the villi and bottom boundary, it avoided a mismatch of boundary conditions and it was observed that very little flow was induced in this region. In multiple-villi simulations, the spacing between the villi was maintained at $0.675L$⁠. The arc length of travel of the tip of each villus was kept constant at $W=Lπ/12$⁠. The geometry is shown schematically in Fig. 3.

On the fluid domain, zero-stress boundary conditions were imposed for the velocity and a zero-normal-gradient condition was imposed for the pressure on all the boundaries. On the villi, a no-slip condition was imposed for the velocity, and a zero-normal-gradient condition was imposed for the pressure. Note that this choice of boundary conditions, along with the interpolation scheme used to impose the immersed boundaries, has been shown to maintain the second-order accuracy of the finite-difference scheme.³⁸.

Although several possible length scales and thus non-dimensionalizations for the Reynolds number exist, we use the definition below to capture effects of regular oscillations and villus width

where ψ is the villus oscillation angular frequency and ν is kinematic viscosity. Note that while the geometry (the shape and spacing of the villi) is not varied, and the oscillation of the villi is sinusoidal, the problem is fully parameterized using Re, the phase angle $ϕ$⁠, and an oscillatory Reynolds number, or the square of the Womersley number $Reosc=ψL2/ν$⁠. These parameters arise using L to nondimensionalize lengths and ψ to nondimensionalize time. We have varied Re by varying the viscosity (equivalent to varying the inverse of the frequency), and as such, both Re and $Reosc$ vary proportionally. We highlight that by presenting our results only in terms of Re, we are not implying that our results are independent of the value of $Reosc$⁠. However, nondimensionalizing in the way chosen and presenting in terms of Re is practical or “natural” from the perspective of applications—the normalized data show the trends that would occur if a given system (biological or synthetic) with fixed geometry was explored as a function of frequency.

In the simulations below, the Reynolds number was varied between $Re=10$ and $Re=6000$ with the following values: 10, 100, 1000, 2000, 4000, and 6000. This was achieved by fixing $ψ=1 rad s−1$⁠, and varying the value of ν (given in units of $L2s−1$⁠).

Three distinct flow regimes were consistently observed under the parameters studied:⁴⁵ 

• Stage 1: Isolated tip vortex enclosing the oscillating villi tip.

• Stage 2: Localized stably oscillating vortex dipole at the same frequency as the villus.

• Stage 3: Vortex shedding and the development of spatiotemporal complexity.

We note that the range of Re of each stage and the value of Re of the transitions between these stages have not been determined exactly; however, the values we report are indicative and clearly show the existence of these three distinct stages. To gain a sense of how the different stages affect flow structure, the generic flow features typical to each stage are shown in Fig. 4; while these results are generated by a single villus, the same structures are observed in multi-villi flow, as shown in Sec. III.

In stage 1, despite the presence of an isolated tip vortex, mixing via symmetric oscillations is limited due to dominating viscous effects that cause the flow to be nearly reversible. In this stage, the flow is quasi-steady; the flow is overwhelmingly defined by the position of the villus, with insignificant effect of the history of the motion, and there is little opportunity for mixing. As inertia increases in stage 2, Fig. 4(b), a portion of angular momentum is preserved from each half-oscillation, resulting in the formation of a pair of counter-rotating vortices, but this dipole remains attached to the tip of the oscillating villus. Further increases in inertia in stage 3, Fig. 4(c), induce vortex shedding away from the tip, resulting in complex nonlinear interactions indicative of mixing.

Since vortex formation is the basic feature of transitional and turbulent flows, knowledge of the common interactions of vortices forms a pathway for understanding the transport and mixing mechanics in flows generated by villus-like structures. Our simulations, described in detail in the following section, show that as inertia increases, oscillatory motions generate sets of vortices that interact in increasingly complex ways. Formation of counter-rotating vortex pairs leads to transport dipoles, which transport fluid over potentially large distances, whereas co-rotating vortices enhance stirring locally through merging. To aid understanding of these behaviors, we give a brief summary of common vortex interactions below; for further details, see, for example, Ref. 46.

By applying conservation of angular momentum, it can be shown⁴⁶ that a pair of inviscid point vortices (dipole), with vorticities $Γ1,2$⁠, located at $x1,2c$⁠, rotate about a fixed point

with constant angular velocity

where b is their constant separation. An equal strength co-rotating vortex pair, $Γ1=Γ2=Γ$⁠, rotates about $Xc$ with angular velocity $Ω=1/(πb2)$⁠. By contrast, an equal strength counter-rotating pair, $Γ1=−Γ2=Γ$⁠, advects along a straight path with speed $U=Γ/(2πb)$⁠. Considering these two cases together, counter-rotating pairs of dissimilar strength orbit along circular paths with large radii and small angular velocity, whereas dissimilar strength co-rotating pairs orbit with large angular velocity and a small mean radius. In all cases, the pair maintain their initial separation distance b from each other.

Since both counter and co-rotating point vortices generate closed streamlines, a vortex pair acts (over short timescales) as a transport mechanism that translates a coherent set of particles. The effects of viscosity ultimately cause co-rotating vortices to undergo an instability called merging that causes the pair to join into a single vortex of larger strength.⁴⁷ This develops via the influence of outer ghost vortices, which draw in diffused vorticity when the pair are sufficiently close. Although the timescale of merging is dependent on the relative strength and distance of the pair,⁴⁸ the merged vortex stirs fluid in a spiral-like pattern, resulting in mixing.⁴⁹ By contrast, counter-rotating dipoles tend to be more stable thanks to their relatively smaller net angular momentum,⁴⁶ facilitating their transport over large distances.

The flow structures and effects described above were found to be present throughout our simulations, and to interact in a variety of complex ways, as we describe in Secs. III–V.

We begin by introducing the effects of transport and mixing generated by the oscillation of a pair of villi-like structures. To explore the effects of interaction between the flows generated by the two villi, we analyze two simple patterns:

1. in-phase motion,

2. anti-phase motion.

These patterns are good test cases for investigating the effects of variable phase because they vary drastically in the way they affect the geometry of the intervillus region, as shown in Fig. 5. Compared to the in-phase pattern, the anti-phase pattern oscillates approximately ten times more fluid from the intervillus region, indicating an increase in inertia from within the intervillus region. This leads to an order more maximum instantaneous mass flux for the anti-phase case, increasing mass flux along the tip line. Therefore, the in-phase and anti-phase patterns serve as useful extreme examples to explore and explain how the dynamic coupling of villi varies with phase difference.

Figures 6 and 7 show snapshots of the vorticity field generated by two oscillating villi as time increases, for a range of different values of Re between $Re=10$ and $Re=6000$⁠, for in-phase and anti-phase motions, respectively. Both in- and anti-phase motions exhibit the three transitional regimes described above.

At low inertial scaling, $Re=10$⁠, the flow exhibits stage 1 behavior and small vortices appear that remain attached to each tip, increasing in intensity and flattening as Re increases.

As the flow transitions to stage 2, by $Re=1000$⁠, an oscillating dipole appears above each villus, for example, the two sample images that are approximately 1 period apart at t = 2.25 and t = 7.25 in Fig. 6 or Fig. 7. Since the sign of the circulation of each vortex is determined by the phase of the villus generating them, the pair of dipoles have the same sign of vorticity for the in-phase pattern and opposite sign for anti-phase pattern. Dipoles are separated by a gap of low vorticity for the in-phase pattern and do not interact. However, the large variation in tip separation of the anti-phase pattern forces the dipole pair to be brought close together, resulting in flattening along the central intervillus axis. Despite this, the vortex dipoles remain fixed to the tip and the vorticity field is mirrored along the intervillus axis throughout the flow in stage 2.

During stage 3 (typified here by the flows at $Re=4000$ and $Re=6000$⁠) of anti-phase forcing, anti-phase dipole flattening results in the collision of detached oppositely signed dipoles. As a result, the vortices lose their partner in favor of the neighboring vortex. This results in two new counter-rotating dipoles of similar strength that have advective components in the vertical direction, as shown in Fig. 8(c). The lower and upper dipoles have downward and upward advective components, respectively. However, this formation occurs while the intervillus volume is increasing, reducing the net velocity of the upper vortex. The flow grows more complex as it develops in time: counter-rotating dipoles of similar strength continue to form and undergo further interactions, leading to transport away from the villi in a variety of directions, as shown in Fig. 8(d). Therefore, the anti-phase pattern displays behaviors that indicate both mixing and transport.

By contrast, in stage 3 with in-phase patterns, vortex shedding begins with the creation of two well-separated dipoles of equal orientation and sign, as shown in Fig. 9(a). The dipole center rotates with high curvature, which causes one of the dipoles to advect toward the outer base while the other rotates into the intervillus region. Since the timescale of advection is similar to the oscillation period, the intervillus dipole interacts with a new dipole being formed by the opposing villus to which it was created. Increasing inertia, or equivalently increasing Re, leads to intensified intervillus–villus dipole interactions that result in the formation of additional counter-rotating dipoles that advect slowly and dissipate quickly. As a consequence, equal strength counter-rotating dipoles are infrequent and rarely advect away from the boundary. Thus, the in-phase pattern displays reduced transport characteristics compared to the anti-phase case.

To compare how the various vortex shedding modes affect transport, we measure the lengths of the particle paths: the total displacement $d(x)$ of a fluid particle (following the flow), over the duration $t=[0,T]$ of the numerical experiment, that started at position $x=x(0)$ at time t = 0,

referred to as a displacement field in what follows. Particles that intersected a villus boundary were removed and had their final position estimated by averaging neighboring particles.

The displacement fields after 10 villus oscillations for in-phase and anti-phase patterns are shown in Figs. 10 and 11, respectively. During Stage 1, both patterns display jet-like behavior along the central intervillus axis that is fed by two entrainment vortices above the tips. Since the anti-phase pattern produces a larger mass flux, increased magnitudes of kinetic energy result in intensified displacements compared with the in-phase case. Since the displacement field is composed of a single jet rather than two, the tip vortices couple to affect displacement over a long timescale.

In stage 2, development of the oscillating vortex pair significantly augments the displacement field for both in- and anti-phase patterns, as shown at $Re=1000$ in Figs. 10(c) and 11(c). Instead of an upward jet of fluid, a small mixing region is identified in the form of a complex displacement field near the tip and shaft regions. Since displacement magnitudes are similar for the anti-phase and in-phase cases, there is an indication of a lack of cross-flow between the villi. Nevertheless, vortex merging causes fluid to mix within the intervillus region, penetrating deeper for the anti-phase case.

In stage 3, the mixing region increases above the tips of the villi, resulting in a mushroom cloud shape, shown in Figs. 10(e) and 11(e). Further increases in Re extend the mixing region in space, both horizontally and vertically, see Figs. 10(f) and 11(f).

Throughout this process, the anti-phase case has consistently larger displacement magnitudes and mixing regions. This supports our proposition that the anti-phase pattern enhances mixing and transport via the collision of opposing dipoles.

We showed above how the flow generated by two oscillating villi exhibits a variety of vortex interactions, dependent on the phase difference of oscillations. Transport was accentuated for the anti-phase pattern due to energy increases associated with intensified intervillus fluxes and formation of additional transport dipoles. To explore how villus coupling extends to larger scales, such as those one might encounter in a biological or artificial gut, in this section, we analyze an array of 15 villi. For ease of comparison with previous results, we generate flow fields for in-phase and anti-phase patterns with the same set of parameter values as before.

Numerically computed vorticity fields are shown in Figs. 12 (Multimedia view), 13 (Multimedia view), and 14 as snapshots in time for a range of different values of Re, with both in- and anti-phase patterns; the results are described in detail below.

In the system with 15 villi, similar instabilities to the two-villi case occur, at Reynolds numbers of the same order of magnitude, as shown in Figs. 12 and 13. Stage 1 behavior consists of singular vortices attached to each tip, which flatten and intensify as inertia increases. During stage 2, in-phase dipoles oscillate independently, whereas anti-phase dipoles flatten on approach, as shown in Fig. 14(a). However, since all but the first and last villi have two neighbors, flattening events occur at twice the rate as the two-villi case, see Fig. 14(b).

However, the increase in dipole-flattening events for the anti-phase pattern causes new emergent behaviors to occur during stage 3. During the first half-oscillation, dipole collisions occur in a similar manner to the two-villi case, resulting in sets of opposing dipole pairs aligned vertically within an opened intervillus region. Over the next half oscillation, the closing of the intervillus region forces the lower dipoles to break apart and flow over each tip, then pairing with a generating vortex. This results in a high strength dipole collision, similar to the two-villi anti-phase case, but at a spacing of two intervillus spaces. The top dipoles of this collision become so strong by $Re=6000$ that vortices advect away from the influence of opening intervillus spaces, resulting in enhanced vertical transport on a short timescale. As the flow grows in time, complex vortex interactions result in increased generation of transport dipoles with inertia.

The in-phase dynamics are comparatively tamer in stage 3, as shown in Fig. 12. Similar to the two-villi case, biased dipoles develop above the tips and rotate sharply into neighboring intervillus spaces. Complex interactions develop only for $Re=6000$ and occur within a smaller vertical region compared to anti-phase oscillations. Formation of similar strength transport dipoles is rare and infrequent, indicating reduced transport on a large timescale. However, the in-phase pattern features regular, synchronized transport of dipoles into neighboring intervillus spaces, indicating possible enhancement of horizontal transport.

Further confirmation of the change in behavior due to dipole pairing and advection can be found in the overall kinetic energy of the flow. Figure 15 shows time histories of the total kinetic energy for the 15-villi cases. Both in-phase and anti-phase cases are plotted. In both phase configurations, the behavior of the total energy is qualitatively and quantitatively similar for $Re=1000$ and $Re=2000$⁠. However, there is a distinct change in its behavior for $Re≥4000$⁠. In the in-phase case, the total kinetic energy oscillates with the villi motion about a mean value that begins to increase from $Re=4000$⁠. In the anti-phase case, the total kinetic energy appears to oscillate about a mean value when $Re≤2000$⁠; however, increasing to $Re=4000$ sees the kinetic energy begin to grow substantially—so much so it indicates these cases may not have saturated or reached statistical stationarity.

We compute displacement fields, defined in Eq. (5), in order to examine the transport properties of the generated flows; the results are shown in Figs. 16 and 17 (showing only the central four villi for clarity). We find that the flow is characterized by a region of high mean horizontal displacement, $d¯(y)$⁠, above the villus tip (⁠$y≈1$⁠), which we refer to as the “peripheral region.”

Qualitatively, stage 1 effects are similar for both in-phase and anti-phase patterns, in which a long-timescale circulation region is formed near each villus tip, as shown in Figs. 16(a) and 16(b) and 17(a), and 17(b). However, since the anti-phase pattern experiences dipole flattening every half period, circulation vortices appear at twice the rate as the in-phase pattern. Both patterns have only limited transport between the base of the villi and the tip in stage 1.

During stage 3, the comparatively large intervillus volume flux associated with the anti-phase pattern greatly enhances the height and displacement magnitudes of the peripheral region, as shown in Figs. 16(e) and 16(f) and 17(e) and 17(f). While the height of the peripheral region for the in-phase pattern is clearly below the villus length (indicated by a sharp decrease in displacement), the region quickly grows in height for the anti-phase pattern, reaching a peak of around 4 villus lengths for $Re=6000$⁠.

To gain a sense of scale, we average the displacement field horizontally across the width of the computational domain

where $X=30L$ for our simulations. The results are shown in Fig. 18.

While the height of the peripheral region is similar for both patterns during stages 1 and 2, there are significant differences during stage 3. Throughout this stage, the in-phase pattern's peripheral region increases in mean displacement with an approximately fixed peripheral region height of around 2 villus lengths, indicating enhanced horizontal transport. Conversely, the anti-phase pattern experiences increases in both mean displacement and peripheral region height, indicating enhanced transport both horizontally and vertically. Only small increases to mean displacement can be observed between the $Re=4000$ and $Re=6000$ cases, indicating a possible limit to the peripheral region height. This apparent saturation of the peripheral region height, which occurs for both in-phase and anti-phase cases further confirms the distinction between stages 2 and 3—only in stage 3 does this height begin to saturate and lose a strong dependence on Re.

Having gained a sense of scale for the peripheral region, we proceed to apply a measure of “mixedness,” to the 15-villi simulations, that permits quantitative comparison of mixing in both the horizontal and vertical directions. A common way to measure mixedness is to measure the distance between an initial seed of tracer particles after allowing the tracers to advect and diffuse for a given time.^50,51 Here, we use the shortest distance between two sets of tracers using a metric proposed in a study of mixing in ciliary locomotion.²⁹ This metric is called the “mixing number,” m(t), and is defined by

where $xi(t)$ and $yj(t)$ are the displacements of the members of the two sets of tracer particles, and N is the number of particles in each set. Smaller m indicates greater mixedness.

In order to compare mixing within the peripheral region across both patterns in the horizontal and vertical directions, we use two tracing schemes. These were formed by halving a 100 × 100 tracer square either vertically or horizontally, forming two stripes, as shown in Figs. 19(a) and 19(b). A summary of these schemes and the directions of mixing they quantify is shown in Table I. In order to compare the speed of mixing in regions with significant transport, the square was aligned to the bottom of the domain and size set to the length of two-villi lengths, the approximate height of the stage 3 in-phase peripheral region.

An example evolution of tracer particles is shown in Fig. 19, for both schemes, for 15 villi in anti-phase motion at $Re=6000$⁠. After two cycles [Figs. 19(c) and 19(d)], we see that while the particles remain largely separated, they begin to overlap in the central regions, decreasing the value of m. Following the fourth cycle [Figs. 19(e) and 19(f)], there is significant overlap, enhanced in part by a transport dipole in the top-right, further decreasing m. By the ninth cycle [Figs. 19(g) and 19(h)], the particles overlap everywhere, indicating a well-mixed solution. Throughout this process, it appears that the vertical scheme [Figs. 19(a), 19(c), 19(e), 19(g)] generates more overlap than the horizontal scheme [Figs. 19(b), 19(d), 19(f), 19(h)], indicating a faster change in the mixing number in the vertical direction, as suggested by all previous analysis.

The mixing number, m(t), was computed for both in-phase and anti-phase patterns, over the full range of Re as a function of time; the results are shown in Fig. 20, with m normalized by m₀, the value of m at the initial particle distribution. From the definition in Eq. (7), it would therefore be expected that m should initially reduce as the two blocks or strips of particles of different color begin to mix. If the mixing process was purely diffusive, it may be expected that $m/m0$ would initially decrease as the two colors mixed, reach a minimum, and then slowly begin to increase again as the particles diffuse over a larger area. The rate of the initial decrease in $m/m0$ would be related to diffusivity of the particles; higher diffusivities resulting in faster mixing would produce a more rapid rate of decrease in $m/m0$⁠. So, cases for which $m/m0$ decreases the most rapidly initially are those which produce the most effective mixing. The oscillatory nature of stage 1 flow causes m(t) to oscillate and decrease slowly, for both the horizontal and vertical schemes in both in- and anti-phase patterns, indicating limited mixing over a timescale of ten oscillations. Although m(t) decreases steadily with time across the entire range of Re for the in-phase pattern [Figs. 20(a) and 20(c)], the rate of mixing sharply increases for the anti-phase pattern during stage 3 [Figs. 20(b) and 20(d)].

By plotting the final mixing number as a function of Re, in Fig. 20(e), we see that the mixing is significantly enhanced by increase in Re, as expected. Furthermore, since the anti-phase pattern has significantly smaller final mixing number than the in-phase pattern, at the same value of Re, for both horizontal and vertical tracer schemes, we can conclude that the anti-phase pattern is the superior pattern for both transport and mixing.

Knowledge of the hydrodynamic response of coupled villi-like actuators across inertial scales has the potential to inspire the design of artificial prototypes for sensing and mixing applications. Previous experimental and numerical studies on villous motility suggest that groups of villi might coordinate to produce efficient peripheral mixing within the transitional flow regime.^25,26 To investigate mixing, we developed a simplified two dimensional computational model consisting of linear arrays of villi-like actuators and used it to explore the flow properties for in-phase and anti-phase oscillations throughout the transitional flow regime. These patterns were chosen because they maximize the difference in intervillus volume change and, thus, formed a good model to explore the effect of a change in coupling.

Simulations involving two villi, Sec. III A, revealed significant differences in vortex dipole formation for the in-phase and anti-phase patterns. The converging tips of the anti-phase pattern produced dipoles of similar strength, which transport fluid away from the villi over large distances for $Re≥2000$⁠. By contrast, in-phase oscillations have a tendency to produce biased dipoles that rotate with high curvature, resulting in comparatively little transport. As a result, displacement fields revealed a peripheral mixing region that grows larger and more complex with increased inertia for the anti-phase pattern.

Augmented hydrodynamic effects were observed in arrays of 15 villi as a result of additional couplings between neighbors. During anti-phase oscillations, the ejection of fluid from an intervillus space coincides with admission of its neighboring spaces, reducing the velocity at which fluid is propelled upward. Since newly formed transport dipoles have a direction of travel away from the periphery, there is a competing effect with local intervillus volume changes. Although these competing effects are approximately equal for $Re≈2000$⁠, dipole advection dominates for $Re>4000$⁠, resulting in the generation of a large peripheral transport layer of height approximately three villus lengths, as indicated in Fig. 18. By contrast, in-phase oscillations rarely produce transport dipoles, resulting in a shorter peripheral region with a height of less than two villus lengths at $Re=6000$⁠, reducing transport.

Mixing timescales were found, in Sec. III B 3, to be consistently faster for the anti-phase pattern in both the horizontal and vertical directions. Since the largest gains in final mixing number for the anti-phase pattern were obtained at $Re≈2000$ (Fig. 20), mixing at high Re is likely driven by the advection of transport dipoles. Thus, despite the competing effects associated with the simultaneous admission and ejection of fluid, the anti-phase pattern was found to be superior to the in-phase pattern for both transport and mixing, across all inertial scales tested.

The results imply a variety of useful control behaviors for mixing and chemical sensing applications. Stage 1 is advantageous for tip sensing as localized circulations can be induced with minimal disturbance to the surrounding fluid. When measurements within the intervillus region are required, stage 2 could be utilized to extend the domain and mix the fluid. Finally, stage 3 could be employed to simultaneously mix, transport, and detect fluid within a larger peripheral region that extends above the villi tips.

Considering robotic applications, we propose that stages could be employed in combination to fulfill more complex goals. For example, an MFC swimming with embedded villi with could program a multi-stage digestion cycle. At the start of the cycle, intake of organic material could be detected using low-power Stage 1 dynamics. Following detection, stage 3 could be employed to boost power output by mixing the fluid evenly near the surface. The system could then progress to stage 2 to balance the needs of circulation, power efficiency, and anti-sedimentation.

In this paper, we investigated transport and mixing in arrays of oscillating villi-like structures in two-dimensional, single-phase transitional flows.

By interpreting the structure of the flow as an ensemble of interacting vortex dipole pairs, we characterized the response to forcing by two villi across a range of length scales, revealing the development of a mixing region from the onset of the turbulent regime. This mixing region was found to significantly increase transport during anti-phase oscillations at high Re through collision of shedding vortex dipoles. With arrays of 15 villi, we found that the anti-phase pattern significantly increases the spatial and temporal scales of the peripheral mixing region across transitional regimes. These increases were driven by the formation of transport dipoles that arise as a result of increased flux from the intervillus regions. Additional couplings of the array of 15 villi resulted in a dramatic growth of the height of the peripheral region, indicating that cross-coupling between villi leads to emergent global behaviors that enhance both transport and mixing.

Our results highlight the effectiveness of coordinated oscillating volumes as a strategy for producing a peripheral mixing region, indicating promise for the development of bioinspired artificial villi. Given that the height of the peripheral region is significantly larger for larger intervillus volume changes, active variable-phase mechanisms are likely to be more advantageous than single-phase strategies. Such gains are likely to be accentuated at higher Re, where complex interactions increase the strength of coupling across networks of villi.

The insight gained from this work may help inform knowledge of digestion in chemical sensors and fluidic mixers. Since the couplings associated with coordinated villus-like motions can be utilized for the control of the domain of fluid circulation, it would seem plausible that these effects would increase the sensitivity of chemical probes. This could be used for the design of future robotics and fluidic-control systems such as MFCs and robotic fish. Of course, further work incorporating facets such as background flows, geometric and kinematic complexity, and variable fluid properties is required to assess the impact of these results in specific applications.

Future work will explore additional patterns, such as groups of villi oscillating in translation as well as (or instead of) in rotation, three-dimensional effects in extended arrays, and soft flexible actuators that more closely represent the villi found in real biological systems. Given the rich variety of dynamics found here, this is expected to yield new hydrodynamic structures that could further enhance the mixing properties of the peripheral region.

The authors have no conflicts to disclose.

Aaron Fishman: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Software (equal); Writing – original draft (equal); Writing – review and editing (equal). Jonathan M. Rossiter: Conceptualization (equal); Supervision (equal); Writing – review and editing (equal). Justin Leontini: Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Software (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review and editing (equal). Martin Homer: Conceptualization (equal); Methodology (equal); Project administration (equal); Supervision (equal); Writing – review and editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.