# Hydrodynamics in a villi-patterned channel due to pendular-wave activity

Rohan Vernekar, Faisal Ahmad, Martin Garic, Dácil Idaira Yánez Martín, Claude Loverdo, Stéphane Tanguy, Clément de Loubens

► To cite this version:

Rohan Vernekar, Faisal Ahmad, Martin Garic, Dácil Idaira Yánez Martín, Claude Loverdo, et al. Hydrodynamics in a villi-patterned channel due to pendular-wave activity. 2025. hal-04988881

HAL Id: hal-04988881

https://hal.science/hal-04988881v1

Preprint submitted on 13 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Hydrodynamics in a villi-patterned channel due to pendular-wave activity

Rohan Vernekar<fnref n="1" />†, Faisal Ahmad¹³, Martin Garic <fnref n="2" />, Dácil Idaira Yánez Martín¹³, Claude Loverdo<fnref n="2" />, Stéphane Tanguy<fnref n="3" /> and Clément de Loubens<fnref n="1" />‡

Inspired by the motility of small intestine, we investigate the flow generated by a propagating pendular-wave along the walls of a symmetric 2D channel lined with elongated villi-like microstructures. The rigid villi follow simple harmonic axial motion, but with a phase lag to their neighbours, resulting in travelling intervillous contractions. We use lattice Boltzmann simulations to resolve the flow around the villi and the luminal channel space, sampling small to moderate Womersley numbers.

We analyse an emergent 'mixing' boundary layer in the lumen, separating two distinct flow regions. A mixing region is formed above the villi, characterised by semi-vortical flow patterns, which travel along with the wave. In the channel centre, unidirectional axial flow emerges, flowing opposite to the imposed wave direction, even in the Stokes flow regime. This behaviour contrasts with canonical wave driven peristaltic flow (Jaffrin & Shapiro 1971). The fluid trapped in the intervillous gaps is forced to take a non-reciprocal path due to the travelling pendular-wave, which breaks time-symmetry in the Stokes regime.

We propose an effective velocity boundary condition, incorporating perturbations due to the presence of the villi. In the inertial regime, the boundary layer shrinks due to dynamic flow confinement arising from increased oscillatory flow inertia. We identify the Stokes to inertial transition, and find phenomenological scaling laws for the boundary layer decay and axial and radial fluid fluxes. Our results not only advance the understanding of transport within the gut but also suggest a novel way for microfluidic flow control in confined channels.

## 1. Introduction

Inspired from motility of smooth muscle tubular structures, such as the gastro-intestinal tract or the ureter, the peristaltic pump stands as one of the most well established bioinspired fluidic system, enabling the efficient pumping of viscous fluids while preventing contamination from downstream to upstream (Esser et al. 2019). The basic mechanism of peristaltic pumping takes advantage of viscous dissipation within the propagating zone of contraction, generating a pressure gradient and subsequently a net flow rate in the direction of the wave. Seminal

<footnote n="1">¹Univ. Grenoble Alpes, CNRS, Grenoble INP, Laboratoire Rhéologie et Procédés (LRP), 38000 Grenoble, France</footnote>

<footnote n="2">²Sorbonne Université, CNRS, Institut de Biologie Paris-Seine (IBPS), Laboratoire Jean Perrin (LJP), F-75005 Paris, France</footnote>

<footnote n="3">³Univ. Grenoble Alpes, CNRS, UMR 5525, VetAgro Sup, Grenoble INP, TIMC, 38000 Grenoble, France</footnote>

<footnote>† Email address for correspondence: rohan.vernekar@orange.fr</footnote>

<footnote>‡ Email address for correspondence: clement.de-loubens@univ-grenoble-alpes.fr</footnote>
theoretical work by Jaffrin and Shapiro (Shapiro *et al.* 1969; Jaffrin & Shapiro 1971), rooted in thin-film flow approximation, has sparked significant interest in extending this concept across various flow regimes, for both biological (Sinnott *et al.* 2017; Amedzrovi Agbesi & Chevalier 2022; Takagi & Balmforth 2011) and engineering applications (Laser & Santiago 2004; Pandey *et al.* 2023).

Beyond peristalsis, biology continues to inspire the development of fluidic systems, particularly in microfluidics. A striking example is the emergence of artificial cilia-based devices that replicate the action of active biological ciliated cells, which measure just a few tens of micrometers (den Toonder *et al.* 2008; Tabata *et al.* 2002; Khaderi *et al.* 2011; Shields *et al.* 2010). These cilia, which line the mucosa of the respiratory system, reproductive tract, and cerebral ventricles (Satir & Christensen 2007; Marshall & Kintner 2008; Shah *et al.* 2009), play a crucial role in clearing viscoelastic mucus (Button *et al.* 2012; Loiseau *et al.* 2020; Choudhury *et al.* 2023). Effective transport in a preferred direction arises from the complex motion of cilia carpets. At the scale of an individual cilium, the beat is inherently non-reciprocal, effectively breaking the time reversibility of low Reynolds number flows. At a larger scales, hydrodynamic interactions promote the emergence of metachronal waves travelling across densely packed cilia carpets (Brennen & Winet 1977; Dauptain *et al.* 2008; Elgeti & Gompper 2013), causing directional fluid pumping. Beyond their role in biological transport, cilia-inspired microfluidic technologies also offer the ability to spatially manipulate flow into distinct regions (Shields *et al.* 2010): within the cilia carpet, vortical flow dominates, enhancing mixing, while above the cilia, long-range unidirectional transport is facilitated by the apparent shear stress generated at the cilia tips.

While cilia-driven transport exemplifies how biological microstructures regulate fluid dynamics at small scales, another striking example of bio-inspiration is found in the small intestines. The inner wall of the small intestine is lined with finger-like or leaf-like (or ridge-like) structures known as villi (figure 1(A)), measuring approximately 500 µm in height. Although their role in passively increasing absorption through surface area augmentation is often overstated (Strocchi & Levitt 1993), in-depth *ex vivo* investigations reveal that these microstructures likely play an active mechanical role in enhancing mixing near the intestinal wall (Westergaard *et al.* 1986; Mailman *et al.* 1990; Levitt *et al.* 1992).

At first glance, one might assume that densely packed intestinal villi function in a manner similar to carpets of ciliated cells in promoting fluid transport. However, from a physiological perspective, their structures are fundamentally different. A cilium is a single cellular extension whose beating pattern is actively controlled by motor proteins, whereas a villus is an elongated multicellular structure composed of an epithelial layer surrounding a network of lymphatic and vascular vessels (Hosoyamada & Sakai 2005). Unlike cilia, which exhibit significant bending flexibility, villi are nearly rigid under physiological flow conditions and do not possess the same degree of freedom for bending motions (Lim *et al.* 2014). Instead, their movement is driven by the contractions of the circular and longitudinal smooth muscles of the intestinal wall (Lentle *et al.* 2013). Numerical simulations have identified two main physical mechanisms by which villi enhance mixing. First, fluid inertia influences the flow field generated by the oscillatory motion of the villi, as the associated Reynolds number (*Re*) is of the order of unity. The periodic forcing induced by these oscillations generates steady streaming flow (SSF) patterns, leading to the formation of a localized micro-mixing layer near the villus tips (Wang & Brasseur 2017; Puthumana Melepattu & de Loubens 2022). This phenomenon arises from time-irreversible convective terms in the Navier-Stokes equations (Riley 2001; Boluriaan & Morris 2003; Wu 2018; Nama 2024; Marmottant 2024), further distinguishing villi-driven mixing from cilia-driven transport. Second, the relative movement between adjacent villi induces periodic variations in the volume of the intervillous gaps. This dynamic expansion and contraction acts as a micro-pumping mechanism, driving
fluid radially in and out, and thereby enhancing radial mixing near the mucosal wall (Lim *et al.* 2015; Lentle *et al.* 2013; Wang & Brasseur 2017; Zhang *et al.* 2020), as illustrated in figure 1(B).

Assuming that villi passively follow the dynamics of mucosal deformations induced by contractions of longitudinal and circular smooth muscle, a variety of contractile motility patterns can be expected (Macagno & Christensen 1980). In this study, we focus on pendular-wave activity, which corresponds to the propagation of longitudinal contractions along the length of the small intestine (Melville *et al.* 1975; Lammers 2005; Lentle *et al.* 2012). To model this motility, we consider villi as rigid microstructures undergoing harmonic oscillations. The underlying wave propagation of the longitudinal muscle layer is manifested by an array of villi oscillating with a constant phase lag $\Delta\phi$ between neighbouring villi. The oscillation angular frequency $\omega$ and number of villi per period ($N = 2\pi/\Delta\phi$) in this array therefore defines the wavelength $L_z$ as well as the wave speed $c$. This framework resembles the metachronal waves observed in ciliary systems (Ding *et al.* 2014; Hall & Clarke 2020). However, unlike cilia-driven flows where viscosity dominates, inertial effects associated with pendular-wave villi contractions are expected to significantly influence the flow-field (de Loubens *et al.* 2013).

For instance, in the rat duodenum filled with ingested water (kinematic viscosity $\nu \approx 7 \times 10^{-7}$ m²/s at body temperature), the Womersley number, defined as $Wo = W/\sqrt{\nu/\omega}$, reaches approximately 0.5 for typical pendular oscillations of ~40 cycles per minute (Lentle *et al.* 2012). Here $W \approx 200$ µm represents the characteristic villus width. In flows driven by oscillatory boundaries, fluid inertia introduces an additional hydrodynamic length scale, the Stokes layer thickness $\delta_{St} = \sqrt{\nu/\omega}$ (Schlichting 1960). This length scale directly influences both the location and strength of the steady streaming flow (SSF), an effect that has important implications for microfluidic mixing (Costalonga *et al.* 2015; Fishman *et al.* 2022). Although SSF is typically associated with inertial effects, it can also arise in Stokes flow conditions when additional degrees of freedom, such as boundary deformability, are present (Marmottant 2024; Cui *et al.* 2024). Cilia-driven flows, for example, may be interpreted as a form of SSF in the Stokes regime, although this classification is not commonly adopted in literature (Riley 2001; Hall & Clarke 2020).

In this study, we demonstrate that Stokes flow reversibility can also be broken by the phase-lagged motion of an array of villi. This mechanism drives irreversible steady streaming flow through the channel, opposite to the direction of the villi-wall pendular-wave. The instantaneous flow field exhibits a boundary layer, the height of which is determined by the nature of the decay of oscillatory flow generated by the villi-wall. Furthermore, we show that as oscillatory inertia increases, the Stokes layer thickness decreases, imposing a dynamic confinement on both the instantaneous flow and the SSF. By mapping the boundary layer height from the Stokes to inertial regime, we identify the critical transition between the two flow regimes. We also derive phenomenological scaling laws to describe the rich flow dynamics observed. We postulate an effective velocity boundary condition for the flow generated just above the villi tips. This model is the first to take into account the perturbations induced by the discrete nature of the villi-patterned wall. We demonstrate the application of this model for long-wavelength travelling contractions in the rat duodenum, the parameters for which are obtained from physiological live-gut motility data.

The remainder of the paper is organized as follows; we define the problem and carry out non-dimensionalization in section 2, and explain the numerical lattice Boltzmann method in section 3. We briefly look at the flow due to non-propagating villi oscillations in section 4, and discuss flow arising from propagating pendular-wave along the villi-wall in section 5. The mixing boundary layer is quantified and analysed in section 6. We then discuss micro-
<figure><img src="image_2.png" /></figure>

(B)

<figure><img src="image_3.png" /></figure>

(C)

<figure><img src="image_4.png" /></figure>

(D)

<figure><img src="image_5.png" /></figure>

Figure 1: (A) Image of leaf-like villi in the duodenal section of the small intestine of rat. (B) The 2D simulation domain consists of *N* oscillating villi with periodic boundaries in the ±z direction, a wall at *r* = 0, and an axis of symmetry at the channel center *r* = *R*. Each villus oscillates around its mean position (orange outline) with an amplitude *a* and a phase lag Δ*φ* relative to its neighbour. The mean positions of adjacent villi are separated by a constant pitch distance *P*. The instantaneous positions (filled green) of the villi show that when one section of the villi-wall contracts, the adjacent section relaxes, thereby pumping and drawing fluid in and out of the intervillous spaces. (C) Illustration of the imposed travelling wave and time-periodic motion of the villi-wall, shown for seven z-axially periodic villi. The filled (green) villi, illustrate the temporal evolution of the boundary at different time fractions *t/T*, going from opaque to transparent. (D) The same seven villi boundary system now seen from the top (z, *t* axes) demonstrating the contractions and expansions of the inter-villus gaps. Both the smallest contraction (filled red rectangle) and the largest expansion (hatched purple rectangle) are seen to travel from right to left (along the −z direction) with time *t*.

mixing and transport, with an emphasis on radial fluid pumping in section 7. We present the rat duodenal flow patterns under physiological conditions in section 8, and present the concluding summary in section 9.

## 2. Problem statement and non-dimensionalization

### 2.1. *Problem statement*

We consider a symmetric 2D channel with diameter $2R$. The channel wall (which we refer to as the villi-wall) is patterned with moving leaf-like or ridge-like villi structures of uniform height $H$ and width $W$, as shown in figure 1. The villi are modelled as rigid structures, and the fluid flow is solved in the intervillous gaps and the channel lumen (Lim *et al.* 2014). The spatial coordinates along the axial and radial dimensions of the channel are given by $\mathbf{x} = (z, r)$, respectively. Along $r$, the channel is bounded by a moving wall at $r = 0$, detailed with $N$ periodic villi, and a symmetry boundary at $r = R$. The domain is periodic along ±z. For the purposes of analysis the domain is parted into two zones; the region patterned with villi where $0 \le r \le H$ is termed the *villi zone*, and the channel space where $H < r \le R$, is termed the *lumen zone*. The villi harmonically oscillate along $z$, in time period $T$, around
their mean axial positions as shown in figure 1(C). The mean villi positions (orange outlines in figure 1) are equispaced by distance *P* along *z*, therefore giving a periodic channel length of $L_z = NP$.

Each individual villus, indexed *i* = 1, 2... *N*, is put into axial translation by propagating contractions of the longitudinal smooth muscle cells of the small intestine (Lentle *et al.* 2012). Consequently, the *i*th villi velocity is modelled as,

$$U_i(t) = \omega a \sin (\omega t + (i - 1)\Delta\phi) \hat{z}, \qquad (2.1)$$

and the *i*th villi axial position is,

$$X_i(t) = \bar{X}_i - a \cos (\omega t + (i - 1)\Delta\phi) \hat{z}, \qquad (2.2)$$

where, $\omega = 2\pi f$ is the circular frequency, $f = 1/T$ is the oscillation frequency, *a* is the oscillation amplitude, *t* is the time, $\Delta\phi$ is the velocity phase lag between adjacent villi, $\bar{X}_i$ is the mean axial villi position and $\hat{z}$ is the unit vector along the axial direction.

Since we model a discrete number periodic villi, we can access only discrete values for the phase lag $\Delta\phi = 2\pi/N$. When $N = 1$, $\Delta\phi = 2\pi$, meaning all villi are moving in-phase. This is equivalent to the case where *N* becomes large ($N \to \infty$), and $\Delta\phi \to 0$, a situation recently studied in 3D by Puthumana Melepattu & de Loubens (2022). When $N = 2$, $\Delta\phi = \pi$, and adjacent villi move with velocities exactly out-of-phase. In the intermediate situation, i.e. $0 < \Delta\phi < \pi$ (or $3 \le N < \infty$), a propagating wave of intervillous contraction is generated along the villi-wall. We term this as the *pendular-wave* activity. This wave travels in the $-\hat{z}$ direction at a speed of $c = L_z/T = 2\pi fP/\Delta\phi$. Figure 1(D) demonstrates how contractions and relaxations along the intervillous spaces (see marked rectangular regions) travel from right to left along *z* with increasing *t*, for a domain with $N = 7$ periodic villi. The bottom walls of the intervillous gaps are given instantaneous velocities that vary linearly between those assumed by their bounding villi.

For forming the channel geometry, we take inspiration from the physiology of the rat duodenum (Lentle *et al.* 2012; Hosoyamada & Sakai 2005; Casselbrant & Helander 2022). Table 1 provides typical values for the gut dimensions and motility in the rat duodenum. From these, we set the geometric ratios $R/H = 5.6$, $H/W = 2.5$, and $P/W = 1.6$, which are held constant for all simulations. Here, the oscillatory amplitude is estimated from the rat duodenal strain rate measurement of $26 \pm 8\% s^{-1}$ during pendular-wave activity (Lentle *et al.* 2012). The fluid is assumed Newtonian, with dynamic viscosity $\mu$ and mass density $\rho$. The flow induced by the villi motion is modelled by the incompressible Navier-Stokes equations given by,

$$\nabla \cdot \mathbf{u} = 0, \qquad (2.3)$$

$$\rho \frac{\partial \mathbf{u}}{\partial t} + \rho \mathbf{u} \cdot \nabla \mathbf{u} = -\nabla p + \mu \nabla^2 \mathbf{u}, \qquad (2.4)$$

where *p* is the pressure and *u* the fluid velocity vector field.

## 2.2. Non-dimensionalization

One observes that there are multiple length and time scales in the problem. In addition to multiple length scales in table 1, we can identify time scales for four physical processes. These are: the villi oscillation timescale $t^\omega = 1/\omega$, the fluid advective timescale $t^{adv} = W/(\omega a)$, the viscous timescale $t^\mu = W^2\rho/\mu$, and time scale imposed by the boundary travelling wave $t^c = W/c = W\Delta\phi/(\omega P)$. In the above we have chosen the villi width *W* as the length scale of the problem, in line with the literature (Puthumana Melepattu & de Loubens 2022; Tatsuno 1973; Kotas *et al.* 2007).

Furthermore, oscillatory boundaries, such as the villi-wall in our case, introduce a
<table><caption>Table 1: Typical values characterizing villi geometry and motility for the duodenum of the rat, from Casselbrant & Helander (2022); Lentle *et al.* (2012); de Loubens *et al.* (2013); Hosoyamada & Sakai (2005).</caption><tr><td>Radius of the small intestine</td><td>R</td><td>~2.8 mm</td></tr><tr><td>Width of the villi</td><td>W</td><td>~200 µm</td></tr><tr><td>Length of the villi</td><td>H</td><td>~500 µm</td></tr><tr><td>Pitch</td><td>P</td><td>~320 µm</td></tr><tr><td>Displacement amplitude</td><td>a</td><td>23 – 51 µm</td></tr><tr><td>Frequency</td><td>f</td><td>0.54 – 0.62 Hz</td></tr><tr><td>Velocity amplitude</td><td>U<sub>0</sub> = ωa</td><td>78 – 197 µm/s</td></tr><tr><td>Fluid density</td><td>ρ</td><td>~994 kg/m<sup>3</sup></td></tr><tr><td>Fluid viscosity</td><td>μ</td><td>0.7 mPa.s</td></tr><tr><td>Stokes layer</td><td>δ<sub>St</sub> = √μ/ωρ</td><td>424 – 454 µm</td></tr><tr><td>Womersley number</td><td>Wo = W/δ<sub>St</sub></td><td>0.44 – 0.47</td></tr></table>

characteristic length scale for viscous dissipation normal to the wall, termed as the Stokes layer, whose length is given by $\delta_{St} = \sqrt{\mu/(\omega\rho)}$ (Schlichting 1960). A priori, it is unclear which length and time scales would dominate the problem at hand. We make the choice for non-dimensionalization as: $\tilde{u} = u/(\omega a)$, $\tilde{p} = pW/(\mu\omega a)$, $\tilde{x} = x/W$ and $\tilde{t} = t\omega$ for the typical velocity, pressure, length and time, respectively, to obtain the non-dimensional momentum equation as,

$$ \frac{1}{\tilde{a}} \frac{\partial \tilde{u}}{\partial \tilde{t}} + \tilde{u} \cdot \nabla \tilde{u} = - \frac{1}{\tilde{a}Wo^2} \nabla \tilde{p} + \frac{1}{\tilde{a}Wo^2} \tilde{\nabla}^2 \tilde{u}. \quad (2.5) $$

Two non-dimensional numbers emerge in the (2.5), namely the Womersley number $Wo = W/\delta_{St}$ (Loudon & Tordesillas 1998) and the reduced oscillation amplitude $\tilde{a} = a/W$. The villi motility in (2.1) can now be similarly re-written in the dimensionless form as,

$$ \tilde{U}_i(\tilde{t}) = \sin(\tilde{t} + (i - 1)\Delta\phi) \hat{z}. \quad (2.6) $$

Three parameters then govern the problem, namely the boundary imposed phase lag $\Delta\phi$, the Womersley number $Wo = W/\delta_{St}$ and the reduced amplitude $\tilde{a} = a/W$, in addition to the constant geometric ratios listed earlier.

We carry out simulations in the limit of small amplitude oscillations, at $\tilde{a} = 0.1$ and $0.2$. The physiological Womersley number in the rat duodenum filled with body temperature (36.5°C) water (table 1) works out to be $Wo \approx 0.46$, and the Stokes layer to be $\delta_{St} \approx 440$ µm. In this study we vary the parameter $\Delta\phi$ from 0 to $\pi$ in discrete steps, and the $Wo$ over three orders of magnitude, going beyond the physiological regime. When $Wo^2 \ll 1$, along with $\tilde{a} \ll 1$, we see that the equation (2.5) reduces to the steady Stokes equation. We therefore expect a Stokes flow regime when $Wo \ll 1$, with the flow governed solely by the travelling wave parameter $\Delta\phi$ introduced through the boundary condition. However, when $Wo \gg 1$ we expect an inertial flow regime, with all terms in (2.5) playing a role.

### 3. Numerical methods

We solve the incompressible Navier-Stokes equations using the lattice Boltzmann method (LBM) in 2D using the D2Q9 lattice, which discretizes coordinate space in two directions and the velocity space in nine directions (He & Luo 1997b; Zhang 2011). We use the two-relaxation-time (TRT) scheme, which significantly reduces the relaxation-time dependent error in the LBM (Ginzburg *et al.* 2008), in order to accurately capture second order flow
phenomena. In the TRT-LBM, “distribution function” (or “population”) g is updated at every lattice node according to the following,

$$g_i^*(t^*, \mathbf{x}) = g_i(t, \mathbf{x}) + \frac{1}{\tau^+} (g_i^{eq+}(t, \mathbf{x}) - g_i^+(t, \mathbf{x})) + \frac{1}{\tau^-} (g_i^{eq-}(t, \mathbf{x}) - g_i^-(t, \mathbf{x})) \quad (3.1)$$

$$g_i(t + \Delta t, \mathbf{x} + \Delta \mathbf{x}) = g_i^*(t^*, \mathbf{x}) \quad (3.2)$$

where, $\mathbf{x}$ is the position coordinate, $t$ the time, and the index $i = 0, 1, \dots, 8$ indicates discrete lattice velocity directions. For the $i^{th}$ direction, $\mathbf{x} + \Delta \mathbf{x} = \mathbf{x} + e_i \Delta t$, where $e_i$ is the lattice velocity vector (Krüger et al. 2017). $\tau^+$ and $\tau^-$ are the two relaxation times associated with symmetric (+) and anti-symmetric (-) parts of the populations which are defined as,

$$g_i^+ = \frac{g_i + g_{\bar{i}}}{2}, \quad g_i^- = \frac{g_i - g_{\bar{i}}}{2}, \quad g_i^{eq+} = \frac{g_i^{eq} + g_{\bar{i}}^{eq}}{2}, \quad g_i^{eq-} = \frac{g_i^{eq} - g_{\bar{i}}^{eq}}{2}, \quad (3.3)$$

where the direction $\bar{i}$ is defined such that $e_{\bar{i}} = -e_i$.

We use the He & Luo (1997a) incompressible equilibrium, which is given below,

$$g_i^{eq}(\mathbf{x}, t) = w_i \rho + w_i \rho_0 \left( \frac{\mathbf{u} \cdot \mathbf{e}_i}{c_s^2} + \frac{(\mathbf{u} \cdot \mathbf{e}_i)^2}{2c_s^4} + \frac{\mathbf{u} \cdot \mathbf{u}}{2c_s^2} \right). \quad (3.4)$$

where $w_0 = 4/9$, $w_{1-4} = 1/9$ and $w_{5-8} = 1/36$ are the lattice direction weights, and $\rho = \rho_0 + \delta\rho$ is the fluid density, taken as the sum of constant and variable parts. For the D2Q9 lattice discretization, $c_s = \Delta x / (\Delta t \sqrt{3})$ is the so called “lattice sound speed”. In simulation (or lattice) unites, we take the nodal distance $\Delta x = 1$, the time step $\Delta t = 1$ and $\rho_0 = 1$. Using the incompressible equilibrium rather than the more popular compressible equilibrium is important to accurately capture the time-integrated steady streaming flow.

The macroscopic variables in the solution are recovered through,

$$\rho = \sum_i g_i \quad \text{and} \quad \mathbf{u} = \frac{1}{\rho_0} \sum_i e_i g_i, \quad (3.5)$$

$$p = c_s^2 \rho \quad \text{and} \quad p_0 = c_s^2 \rho_0, \quad (3.6)$$

where $p$ gives the fluid pressure, and $p_0$ is the datum pressure. Through Chapman-Enskog analysis, the kinematic viscosity is related to the relaxation time as $\nu = c_s^2(\tau^+ - \Delta t/2)$ (Ginzburg et al. 2008).

The LBM algorithm follows two simple steps; first, “collision” where (3.1) is evaluated at some intermediate time $t^* > t$, followed by the “streaming” equation (3.2), where the post-collision populations advance to their neighbouring nodes at the end of $t+\Delta t$. This completes one time step $\Delta t$ at the end of which the macroscopic variables $\mathbf{u}$ and $p$ are computed from (3.5) and (3.6). These are then used to compute equilibrium distribution $f^{eq}$ from (3.4) for the next time step.

### 3.1. Moving boundary condition

The complex boundary in this study demands a higher-order LB scheme for treating moving boundaries in order to obtain an accurate solution for the instantaneous and steady streaming flows around the villi and in the lumen. It is equally important to adopt a robust “fresh” node treatment alongside the scheme, in order to ensure a stable solution, and prevent degradation of accuracy. Fresh nodes are those lattice nodes that transition from inside the villi into the fluid zone as the villi boundary moves, and where $\rho$ and $\mathbf{u}$ need to be guessed. We therefore use a second-order accurate interpolated bounce-back scheme (IBB), alongside an iterative procedure to re-fill values for fresh fluid nodes (Ginzburg et al. 2008; Chen et al. 2014).
The IBB scheme is adapted from the family of linear interpolation (LI) schemes that are computationally local at a lattice node (Ginzburg et al. 2023). If the boundary node on the fluid side is at $\mathbf{x}_F$, and the villi wall lies at $\mathbf{x}_W = \mathbf{x}_F + q\mathbf{e}_i$, where $q$ is the fractional distance to the wall along direction $i$, the unknown population along $\bar{i}$ to be streamed from the wall to node at $\mathbf{x}_F$ is computed as,

$$g_{\bar{i}}(\mathbf{x}_F, t + \Delta t) = a_1 g_i^*(\mathbf{x}_F, t^*) + a_2 g_{\bar{i}}^*(\mathbf{x}_F, t^*) + \mathcal{A}_i + \mathcal{W}_i, \quad (3.7)$$

$$\mathcal{A}_i = a_3 \left( I^{FF} g_i(\mathbf{x}_F, t + \Delta t) + (1 - I^{FF}) g_i(\mathbf{x}_F, t) \right), \quad (3.8)$$

$$\mathcal{W}_i = a_0 \frac{w_i}{c_s^2} \rho_0 \frac{\mathbf{u}_W(t) + \mathbf{u}_W(t + \Delta t)}{2} \cdot \mathbf{e}_i, \quad (3.9)$$

where $a_{0-3}$ are interpolation coefficients and $\mathbf{u}_W$ is the wall velocity. We set $I^{FF} = 1$ when $\mathbf{x}_{FF} = \mathbf{x}_F - \mathbf{e}_i \Delta t$ is a fluid node, whereas $I^{FF} = 0$ when $\mathbf{x}_{FF}$ is a solid node (i.e. lies within the villi).

The interpolation weights are set according to the following,

$$a_0 = \frac{3}{5} \left( \frac{4}{1 + 2q} \right), \quad a_1 = a_0 \left( \frac{1}{2} + q \right) - 1, \quad a_2 = 1 - \frac{1}{2}a_0, \quad a_3 = 1 - (a_1 + a_2). \quad (3.10)$$

Note that in (3.10), the prefactor for $a_0$ is set to $(3/5)$ for reasons of simulation stability and to minimising velocity oscillations due to fresh node transitions (though this prefactor value can be varied between [0, 1]). For all reported simulations (except for those presented section 8) we take 20 lattice nodes for the width $W$ of each villi and 32 lattice nodes for the pitch $P$.

### 3.1.1. Fresh node treatment

The process of re-filling all $g_i$ values at the nodes that are uncovered into the fluid side (as the villi move) is called as fresh node treatment. This is a major source of error as well as of numerical velocity oscillations in our computations, and therefore needs to be handled in a robust manner (Chen et al. 2014; Ginzburg 2025). We adopt the local iteration refill (LIR) procedure, with slight modifications, which is highly effective at countering these shortcomings (Tao et al. 2016). In the LIR we first identify groups of connected fresh nodes (at least one fresh node neighbour along $\mathbf{e}_i$), and carry out the following procedure:

(i) Execute local collision step for all link-wise fluid neighbour nodes identified for a fresh node group using (3.1).

(ii) Fresh nodes are partially re-filled by streaming-in post-collision populations from their neighbour nodes using (3.2) as $g_i^*$.

(iii) Fill-in the remaining unknown populations $g_i^*$ using interpolated bounce-back from (3.7).

(iv) Compute temporary populations moments $(\rho^*, \mathbf{u}^*)$ from (3.5).

(v) Then compute $g^{eq,*}(\rho^*, \mathbf{u}^*)$ at the fresh nodes from (3.4).

(vi) Carry out local collision step for all fresh nodes using (3.1).

Steps (ii)-(vi) are repeated for 5 inner-iterations, which has been shown to give sufficiently converged values of $(\rho^*, \mathbf{u}^*)$ (Marson et al. 2021). During the inner-iterations, we update $\rho, \mathbf{u}$ and populations $g^{eq}, g$ only at the fresh nodes (and not at their link-wise neighbours).

### 3.2. Numerical convergence

We start computations with the fluid in the channel at rest, and with advancing simulation time, obtain continuous evolution of the instantaneous and time-integrated steady streaming flow (SSF). At the end of every oscillation period $T$, we compute the developing SSF $\mathbf{u}^{ss,*}$.
by numerically integrating the LBM solution over all time steps of the last competed time period. At the end of a time period, say when $t = T'$, this can be written as,

$$ \mathbf{u}^{ss, \bullet}(\mathbf{x})|_{T'} = \frac{\Delta t}{T} \sum_{i}^{T/\Delta t} \mathbf{u}(\mathbf{x}, T' - T + i\Delta t). \quad (3.11) $$

This developing SSF velocity field is then compared with the developing SSF field computed at the end of the previous time period, say when $T'' = T' - T$, using the L2 norm as,

$$ E_{L2} = \sqrt{\frac{\sum_x (\mathbf{u}^{ss, \bullet}(\mathbf{x})|_{T''} - \mathbf{u}^{ss, \bullet}(\mathbf{x})|_{T'})^2}{\sum_x (\mathbf{u}^{ss, \bullet}(\mathbf{x})|_{T''})^2}}. \quad (3.12) $$

When $E_{L2}$ falls below 1%, we consider that the SSF has converged, and the simulation is set for termination after running for one additional time period of $t = T' + T$. In this converged time period, the instantaneous velocity field, $\mathbf{u}(\mathbf{x}, t)$, and the SSF, $\mathbf{u}^{ss}(\mathbf{x})$, are saved for the given simulation parameters.

## 4. Non-propagating oscillations

As stated in section 2, we carry out simulations by varying the non-dimensional parameters namely the phase lag $\Delta\phi$, and the Womersley number $Wo$, for two values of reduced oscillation amplitude $\tilde{a} = 0.1$ and $0.2$, while the remaining geometric ratios are held constant. Change in $\Delta\phi = 2\pi/N$ results in an associated change in the periodic axial domain length $L_z = NP$, and a change in $Wo$ can be viewed as a change in the oscillation frequency of the villi. In this section we briefly tackle the limiting cases $\Delta\phi = 0$ and $\pi$, which result in non-propagating contractions. The flow velocity contours and streamlines presented in the remainder of the paper are shown for $\tilde{a} = 0.2$, and are rescaled by the characteristic villi oscillation velocity $U_0 = \omega a$, unless otherwise stated.

### 4.1. Instantaneous flow

When $\Delta\phi = 0 = 2\pi$ or $\Delta\phi = \pi$, there is no propagating wave of intervillous contractions along the villi-wall. The figure 2 shows instantaneous flow-fields for $\Delta\phi = 0$, at three time instances, for (A) $Wo = 0.5$ and (B) $Wo = 2.82$. This is the 2D equivalent of the synchronised oscillations of an infinite array of villi recently studied by Puthumana Melepattu & de Loubens (2022). Here the flow (both in the villi and lumen zones) simply follows the villi-wall oscillations. At increased $Wo$, the velocity field becomes localized near the villi, indicating dynamic radial confinement of the flow due to increased oscillatory inertia in the fluid and a decreasing Stokes layer $\delta_{St}$. Here, due to increased inertia the flow in the lumen zone above the villi tips begins to lag from that of the oscillating villi. The opposing axial flow in lumen relative to that within the villi zone leads to momentary re-circulatory flow patterns just above the villi tips, as seen in the second and third panels of figure 2(B). When $\Delta\phi = \pi$, neighbouring villi oscillate exactly out of phase, as shown for three time instances in figure 3, generating counter-rotating vortical flow in the domain, for both (A) small $Wo = 0.16$ and (B) larger $Wo = 2.82$. With increased $Wo$, these vortical flow structures become localized near the villi tips (similarly to the case $\Delta\phi = 0$) due to radial confinement by the Stokes layer.

### 4.2. Steady streaming flow

We compute the steady streaming flow-field (SSF) $\mathbf{u}^{ss}$ as the time-integrated mean flow, using (3.11). These are plotted in figure 4, for the non-propagating phase lags (row-wise) (A) $\Delta\phi = 0$ and (B) $\Delta\phi = \pi$. For a given $\Delta\phi$, we show panels for three increasing $Wo = 0.16, 1.58$ and $5.0$.
10

<figure><img src="image_6.png" /><figcaption>Figure 2: Instantaneous flow-field for the $\Delta\phi = 0$ (synchronous villi oscillations) and $\tilde{a} = 0.2$, with no intervillus contractions, for two different Womersley numbers (row-wise), (A) $Wo = 0.5$ and (B) $Wo = 2.82$, for three time-fractions (column-wise) in one period $T$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 3: Instantaneous flow-field for the case $\Delta\phi = \pi$ and $\tilde{a} = 0.2$, resulting in maximum intervillous contraction and expansion, for two different Womersley numbers (row-wise), (A) $Wo = 0.5$ and (B) $Wo = 2.82$, for three time-fractions (column-wise) in one period $T$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Figure 4: Steady streaming flow-fields ($u^{ss}$) for $\tilde{a} = 0.2$, for the cases with (A) $\Delta\phi = 0$ and (B) $\Delta\phi = \pi$. The panels along the columns are at increasing $Wo = 0.16, 1.58$ and $5.0$. The colour field maps the magnitude of non-dimensional steady streaming velocity ($|u^{ss}|/U_0$).</figcaption></figure>

(column-wise). The SSF is characterised by counter-rotating vortices just above the villi tips. The strength of these vortices is seen to be at least an order of magnitude smaller than villi characteristic velocity $U_0$, except at the highest $Wo = 5.0$. For the synchronised oscillation case $\Delta\phi = 0$, an outer region of weak unidirectional steady flow develops in the lumen with increasing $Wo$, in agreement with the Rayleigh streaming theory (Puthumana Melepattu & de Loubens 2022; Tatsuno 1973). A decreasing $\delta_{St}$ also causes a localization of the flow in this outer streaming region nearer the villi zone with increasing $Wo$. No such outer streaming region develops for $\Delta\phi = \pi$, though the extent of the vortical flow reaches further into the lumen with increasing $Wo$. The key point here is that for non-propagating contractions, $\Delta\phi = 0$ and $\pi$, irreversible SSF manifests as an inertial phenomenon, and the $u^{ss}(x) \approx 0$ in the Stokes flow regime when $Wo \ll 1$.

## 5. Propagating pendular-wave

In this section we deal with the flow dynamics generated by the propagating contractions (and relaxations) of the intervillous spaces along $-\hat{z}$, as illustrated in figure 1(D). These involve cases simulated when $0 < \Delta\phi < \pi$, for increasing $Wo$, with $\tilde{a} = 0.1$ and $0.2$.
<figure><img src="image_9.png" /><figcaption>Figure 5: Snapshots of the instantaneous flow field for $\bar{a} = 0.2$ and $Wo = 0.16$, at two time-fractions (column-wise) in one period $T$, for (A) $\Delta\phi = \pi/2$ and (B) $\Delta\phi = \pi/4$. The dashed (magenta) line marks the approximate separation between the mixing layer and the advecting layer. Note the absence of the advected layer in (B). See supplementary movies 1 and 2.</figcaption></figure>

### 5.1. *Instantaneous flow*

We show two typical examples of the instantaneous flow-fields for low $Wo = 0.16$ in figure 5 and for moderate $Wo = 2.82$ in figure 6. Each individual figure, shows the two time instances (column-wise) of the flow-field, for two separate $\Delta\phi$ values (row-wise). We assert that for $Wo = 0.16$, the flow physics falls well within the Stokes flow regime, while that for $Wo = 2.82$ is in the inertial regime.

From the two figures, we can immediately identify some striking common features of the flow-field structure. We notice the emergence of a boundary layer at distance above the villi tips (approximately marked by dashed magenta line), that separates the flow into two distinct regions. The velocity field between the villi (grey) and this boundary is characterised by asymmetric counter-rotating vortical flows of unequal strength. These vortical flows originate in the contracting intervillous gaps, and terminate in those that are expanding. We call this region the *mixing layer* and the associated boundary thickness as the *mixing boundary layer (MBL) height*, $\ell$. Above the MBL height, until the radial symmetry boundary at $r = R$ we observe near-uniform unidirectional flow moving along $+\hat{z}$. We term this region as the *advecting layer*.

Curiously, the vortical structures in the mixing layer are seen to propagate along $-\hat{z}$, carried by the imposed propagating pendular-wave along the wall, while the flow in the advecting layer goes opposite, along $+\hat{z}$. For the case $\Delta\phi = \pi/4$ and $Wo = 0.16$ (figure 5(B)), we see that the mixing layer extends all the way up to $r = R$, and the advecting layer disappears. In general, the height of the MBL increases with decreasing $\Delta\phi$, until it reaches the channel center-line. Additionally, for a given value of $\Delta\phi$, the MBL height shrinks with an increase
<figure><img src="image_10.png" /><figcaption>Figure 6: Snapshots of the instantaneous flow-field for $\bar{a} = 0.2$ and $Wo = 2.82$, at two time-fractions (column-wise) in one period $T$, for (A) $\Delta\phi = \pi/2$ and (B) $\Delta\phi = \pi/4$. The dashed (magenta) line marks the approximate separation between the mixing layer and the advecting layer. Note that the mixing layer height increases with decreasing $\Delta\phi$. And for a given $\Delta\phi$, it is seen to decrease with increasing $Wo$ (compare with Figure 5). See supplementary movies 3 and 4.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 7: Evolution of the radial velocity component $u_r/U_0$, for $\Delta\phi = \pi/2$ and $\bar{a} = 0.2$, measured at height $r/H = 2$. It is plotted for increasing time-fractions $0 \le t/T < 1$ (opaque to transparent) over the axial domain $z/L_z$, for (A) $Wo = 0.50$ and (B) $Wo = 2.82$. We see travelling velocity wave kinematics, moving from right to left ($-z$ direction) in both panels.</figcaption></figure>

in $Wo$, as can be observed when comparing figures 6 with 5. The decrease of $\ell$ with $Wo$ is shown in supplementary movie 5. For example, in the aforementioned case (at $\Delta\phi = \pi/4$) we see that the MBL shifts from the center-line for $Wo = 0.16$ to approximately halfway between the center-line and villi tips for $Wo = 2.82$. This decrease in MBL is a consequence of the decreasing $\delta_{St}$ with increasing $Wo$, which dynamically confines the flow towards the oscillating villi-wall boundary (Schlichting 1960), as will be discussed later in section 6.

The observed flow dynamics exhibits an intriguing phenomenon; an advecting layer that
axially pumps fluid (along +ẑ) counter to the imposed wall-propagating wave, supported by an underlying mixing layer which travels along with the imposed wave (along −ẑ). Notably, this behaviour persists even in the Stokes flow regime (*Wo* ≪ 1, figure 5), demonstrating symmetry-breaking by the propagating wave. The effect of propagating wave along the villi-wall is felt throughout the lumen zone, irrespective of the MBL. Figure 7 shows the measured radial velocity *u*<sub>r</sub> at *r* = 2H in the lumen, for (A) *Wo* = 0.5 and (B) *Wo* = 2.82, at different time-fractions *t*/*T*. Right-to-left travelling wave kinematics is clearly seen in these plots, with a decrease in wave amplitude with increased *Wo*. A similar behaviour, but with a non-zero positive mean, is also seen for the axial velocity component *u*<sub>z</sub> (not shown here). This behaviour of the instantaneous velocity is seen both in the mixing and the advecting layer, and the wave amplitude is seen to monotonically decreases with increasing radial distance *r*.

5.2. Steady streaming flow

The axial fluid pumping effect is not limited to within the advecting layer, but occurs
throughout the lumen including in the mixing layer. This is demonstrated by the figure 8
with (row-wise) subplots (A), (B) and (C) showing the time-integrated SSF computed for
increasing Wo = 0.16, 1.58 and 5.0, respectively. For each Wo we show five panels for
decreasing values of Δφ (column-wise). The flow is shown over two neighbour villi (at
their axial mean positions) in each simulated system, where the color-field plots the axial
component (u_z^{ss}) of the SSF velocity. We see that the SSF flow is axially periodic about the
intervillous distance (P), around a single villus, in all the panels. The dashed line in the plots
shows the approximate height of the MBL seen in the instantaneous flow-field for each case.

For Wo = 0.16 and 1.58, flow in the lumen is highly uniform in the +r direction, except in and around villi zone. Nearer the villi, flow becomes serpentine (along +ẑ), descending down through the intervillous spaces and curving above the villi tops. We see a small SSF vortex just above the intervillous space, as this serpentine flow meets the uniform flow above. With an increase in Wo = 5.0, the flow is no longer radially uniform in the lumen and its strength gets greatly confined nearer the villi tips, with near zero net flow near the channel center-line. This effect is again a consequence of the increased oscillatory inertia, and the emergence of a flow confining Stokes layer thickness δ<sub>St</sub>. Also, the irreversible SSF flow through the lumen remains a function of Δφ, at all Wo.

We note that the irreversible flow patterns seen in panels of figure 8 cannot be directly
correlated with the emergence of mixing layer boundary discussed in section 5.1. In the
Stokes flow regime (Wo = 0.16), the SSF flow maintains equal strength in both the mixing
and advecting layers for a given (Δφ, til). In the inertial flow regime (Wo = 5.0), the strength
of the SSF decreases radially. This radial decrease in SSF is not dependent on the location
of mixing to advecting layer transition. That is to say, the change in the MBL height ℓ with
Δφ for any Wo, does not appear to leave any signature on the radial variation of the SSF
magnitude (as evidenced by the dashed lines in figure 8).

5.3. Stokes flow regime

The existence of irreversible fluid motion in the Stokes flow regime (Wo ≪ 1) under simple harmonic oscillations of a rigid boundary presents a fundamental departure from the expected reversible behaviour of viscosity-dominated flows. Counter-intuitive to momentum transport *via* travelling waves, and contrasting with the canonical peristaltic flow theory (Jaffrin & Shapiro 1971), the direction of irreversible pumping here is opposite to that of the imposed wall propagating wave. This phenomenon merits separate discussion in the present section.

At first glance, the irreversible axial flow pumping by the villi-wall studied here is
reminiscent of that generated by a propagating waves along dense oscillating ciliary arrays
<figure><img src="image_12.png" /><figcaption>Figure 8: Steady streaming flow-field ($u^{ss}$) streamlines for $\tilde{a} = 0.2$, plotted around a pair of adjacent villi for three increasing Womersley numbers, (A) $Wo = 0.16$, (B) $Wo = 1.58$, and (C) $Wo = 5.0$ (row-wise). For each $Wo$, panels for five decreasing $\Delta\phi$ values (column-wise) are shown. The colour-map plots the axial component of the steady streaming flow (SSF), $u_z^{ss}$. Note that the SSF pattern shows axial periodicity over the intervillous distance $P$, for all cases. The dashed (orange) lines in each panel show the approximate mixing layer height $\ell$ seen in the respective instantaneous flow-fields.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 9: Plots for the (A) Axial oscillatory, $Q_z^{\text{os},v}$, and (B) axial steady-streaming, $Q_z^{\text{ss},v}$, fluxes measured within the villi zone ($0 \le r/H \le 1$), plotted against the phase lag $\Delta\phi$, for increasing $W_o$. Filled markers (coloured) represent cases for $\tilde{a} = 0.2$ and empty markers (black) for $\tilde{a} = 0.1$. Appropriate rescaling causes a reasonable data collapse onto master curves in both plots. In (A), the dashed (blue) curve is obtained from least-squares fitting of the data, with the fit parameters $A = 0.024$ and $B = 0.287$. In (B), the solid (green) curve is derived from geometric considerations in section A.1. Both the curves are functions solely of $\Delta\phi$, demonstrating an axial flux dependence in the villi zone arising purely from geometric motion due to the pendular-wave.</figcaption></figure>

(Ishikawa 2024; Khaderi et al. 2011; Ding et al. 2014; Hall & Clarke 2020). Specifically, those undergoing antiplectic metachronal waves, where the direction of wave propagation is opposite that of fluid pumping by the cilia array. Shields et al. (2010) and Ding et al. (2014) also find that the propagating motility of the cilia arrays divides the flow into a mixing layer, which remains limited to below the cilia tips, and a fluid pumping layer above the cilia. However, the individual cilia are deformable and trace non-reciprocal effective and recovery strokes in their oscillation, which breaks the time-symmetry of the Stokes equation and generates irreversible flow or motion. This is clearly not true for the villi, which are rigid (de Loubens et al. 2014) and trace reversible harmonic oscillations.

The non-reciprocal motion of the villi-wall is therefore more subtle. Here, we shall elucidate the cause of the symmetry breaking in Stokes flow, which comes from the travelling contractile wave along its axis, leading to irreversible SSF. In this paper we shall treat all simulated cases at $W_o = 0.16$ as within the Stokes flow regime. The flow solutions at $W_o = 0.16$ show negligible difference from those at $W_o = 0.28$. In order to verify that these simulations indeed fall within the Stokes regime, we carry out additional simulations with the villi-wall boundary (and its instantaneous velocity condition) held stationary, at given selected time-steps. The converged flow from these simulations is compared with the instantaneous flow solution at $W_o = 0.16$, at matching time-steps, and shows negligible differences in the velocity field.

We recall the *Scallop Theorem* articulated by Purcell (1977), that the back and forth reciprocal motion of the shell of a single-hinged scallop will never generate net irreversible motion or flow, in the Stokes flow regime. However, this reciprocity could be circumvented by introducing an additional degree of freedom to the motion, such as allowing the scallop to periodically rotate about its hinge alongside the opening and closing of its shell. This combined motion will enable the scallop to swim in the Stokes flow regime (Lauga 2011). Using this idea, Najafi (2004) proposed a simple three sphere 1D swimmer, where the gaps connecting the spheres undergo cyclic non-reciprocal contraction-relaxations, while the
<figure><img src="image_14.png" /><figcaption>Figure 10: Axial velocity $u_z/U_0$ measured at the villi tips, $r/H \approx 1$, from simulations at $Wo = 0.16$ and $\tilde{a} = 0.1$, for (A) $\Delta\phi = \pi/2$ and (B) $\Delta\phi = \pi/4$ (dot-dashed, blue) compared with corresponding effective velocity model $\tilde{u}_z^e$ from equation 5.3 (solid, red), for four increasing time-fractions.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 11: Radial velocity $u_r/U_0$ measured at the villi tips, $r/H \approx 1$, from simulations at $Wo = 0.16$ and $\tilde{a} = 0.1$, for (A) $\Delta\phi = \pi/2$ and (B) $\Delta\phi = \pi/4$ (dot-dashed, blue) compared with corresponding effective velocity model $\tilde{u}_r^e$ from equation 5.4 (solid, red), for four increasing time-fractions.</figcaption></figure>

spheres trace reciprocal motion. Revisiting the figure 1(D), it is clear that the same non-reciprocal sequence is carried out by any two adjacent intervillous gaps (and their bounding villi), as by those in the three sphere simple swimmer (Najafi 2004; Najafi & Golestanian 2005). The villi-wall system in the present study can be viewed as an array of such non-reciprocally contracting (and relaxing) gaps, where the additional degree of freedom is enabled by the phase lag $\Delta\phi$. Referring to figure 1(D), the phase-lagged motion simply means that any downstream (along +z) intervillous gap always reaches its minimum (or maximum) size ahead of its upstream neighbour, and this action is never reversed, resulting in non-reciprocal motion.

The non-reciprocal motion of adjacent intervillous gaps transfers fluid from an upstream
(along $-z$) gap to its downstream neighbour. The villi-wall boundary system therefore irreversibly pumps fluid, in the axial direction through the villi zone ($0 \le r \le H$). From the simulations, we can measure the irreversible axial flux within the villi zone using numerical integration as,

$$Q_z^{ss,v} = \frac{1}{L_z} \int_0^H \int_0^{L_z} u_z^{ss}(z,r) dz dr \qquad (5.1)$$

and the oscillatory flux as,

$$Q_z^{os,v} = \frac{1}{2TL_z} \int_0^H \int_0^{L_z} \int_t^{t+T} |u_z(t,z,r) - u_z^{ss}(z,r)| dtdzdr. \qquad (5.2)$$

Figure 9 shows the measured (A) oscillatory axial ($Q_z^{os,v}$) and (B) steady streaming axial ($Q_z^{ss,v}$) fluxes within the villi zone, for various $Wo$ (not limited to the Stokes regime). The oscillatory flux in 9(A) shows a reasonable collapse of the form $\sim \cos(\Delta\phi)$, when rescaled by the reference flux $U_0H$. Here, we see that the data for $\tilde{a} = 0.1$ shows functional collapse at a slightly different level than that for and $\tilde{a} = 0.2$. We fit the equation $A \cos(\Delta\phi) + B$ for all the data in the plot, as a definitive geometric dependence for the oscillatory flux $Q_z^{os,v}$ arising from the imposed motion of the villi-wall has yet to be identified. The least-square fitting coefficients are found to be $A = 0.024$ and $B = 0.287$ used to plot the dashed (blue) curve shown in figure 9(A).

For the irreversible SSF flux in figure 9(B), we see a remarkable data collapse when the flux is rescaled by $\tilde{a}U_0H$. Here, the data begins to show deviations due to increasing inertial effects only for high $Wo = 5.0$ for $\tilde{a} = 0.1$. Taking into account solely the geometric time-averaged motion of the fluid trapped in intervillous gaps, we obtain the scaling for $Q_z^{ss,v} \approx (\omega a^2 H \sin(\Delta\phi)) / (2P)$. The details for the scaling are provided in appendix A.1. This scaling equation gives an excellent description of the collapsed data, as shown by the solid (green) curve in figure 9(B).

The measured data for both $Q_z^{os/ss,v}$ shows minimal dependence on $Wo$, indicating that the flow in the villi zone remains well-approximated by the Stokes regime, even at higher $Wo$. This behaviour arises from the significant geometric confinement of the flow within the intervillous gaps, characterized by a small $(P-W)/H$. As we shall see in section 7, this is not the case in the lumen zone.

### 5.3.1. Effective velocity boundary condition

In order to enable an analytical solution or simulation of the flow in the lumen, without explicitly accounting for the villi, we explore an effective velocity boundary condition (at the villi tips) for the lumen zone. Towards this end, we carry out a Fourier decomposition of the instantaneous velocity signal measured at $u(t,r \approx H)$ for $Wo = 0.16$, just above the villi tips. From the analysis, we find that first order modes of the velocity appear at length scales of the axial domain $L_z$, and multiple second order modes appear at length scales of the villi $P$. We therefore propose the following component-wise non-dimensional effective velocity
<figure><img src="image_16.png" /><figcaption>Figure 12: Coefficients for effective velocity boundary condition in (5.3) and (5.4), as functions of $\Delta\phi$, in the Stokes limit ($W_0 = 0.16$). Empty and filled markers are for $\tilde{a} = 0.1$ and $0.2$, respectively. $c_0$ and $d_0$ are obtained as averages of numerical integration over $z$ and $t$ of the axial and radial velocities above villi tips ($r = 1.02H$). The remaining coefficients are obtained by averaging the least-square fits of the model to the velocity measured above villi tips from the simulations at different times-steps. The standard deviations for the fits are smaller than the size of the markers, and hence not shown. Heuristic rescaling of the coefficient values is carried out to obtain reasonable fitting estimates in the physiologically relevant limit $\Delta\phi \to 0$.</figcaption></figure>

boundary model, $\tilde{u}^e$, for the system,

$$ \tilde{u}_z^e(t, r \approx H) = c_0 + c_1 \sin \left( 2\pi \left( \frac{t}{T} + \frac{z}{L_z} \right) - \frac{\Delta\phi}{2} \right) - c_2 \sin \left( 2\pi \frac{z}{P} \right) \cos \left( 2\pi \left( \frac{t}{T} + \frac{z}{L_z} \right) - \frac{\Delta\phi}{2} \right), \quad (5.3) $$

$$ \tilde{u}_r^e(t, r \approx H) = d_0 - d_1 \cos \left( 2\pi \left( \frac{t}{T} + \frac{z}{L_z} \right) - \frac{\Delta\phi}{2} \right) - d_2 \cos \left( 2\pi \frac{z}{P} \right) \cos \left( 2\pi \left( \frac{t}{T} + \frac{z}{L_z} \right) - \frac{\Delta\phi}{2} \right), \quad (5.4) $$

where $c_{0-2}$ and $d_{0-2}$, are model fitting coefficients that will depend on $\Delta\phi$ and $\tilde{a}$, as well on the geometric ratios of the domain. The functional forms of (5.3) and (5.4) are plotted in figures 10 and 11, respectively, and show an excellent match with the velocity components measured from simulations.

The coefficients $c_0$ and $d_0$ are obtained by taking the time and spatial average of the respective velocity signals from the simulations, for different $\Delta\phi$ and $\tilde{a}$. The remaining fit coefficients ($c_{1,2}$ and $d_{1,2}$) are obtained by least-squares fitting of the model equations to the velocity signal, for all time values, and subsequently taking their mean. These fit coefficients are shown in figure 12 as functions of $\Delta\phi$ for $W_0 = 0.16$. The fit coefficients show complex dependence on $\Delta\phi$ and $\tilde{a}$ in the Stokes flow regime. As the dependence on $\Delta\phi$ and $\tilde{a}$ is not
explicitly clear, we perform a heuristic scaling of the parameters to achieve a reasonable approximation in the limit $\Delta\phi \to 0$. The limiting values are marked (in red) on the y-axis in all panels in figure 12. The coefficients obtained are, $c_0 = 0.28\tilde{a} \sin(\Delta\phi)$, $c_1 = 1.12 - 0.5\Delta\phi$, $c_2 = 0.83\tilde{a}^{-1/4} \tan^{-1}(\Delta\phi/2)$, $d_0 = 0$, $d_1 = 3.23 \tan^{-1}(\Delta\phi/2)$ and $d_2 = 1.84\tilde{a}^{-1/4} \sin(\Delta\phi/2)$. We see that $c_0$ is non-zero, which is necessary for generating irreversible axial flow in the Stokes regime, similar to the effective velocity model for antiplectic ciliary waves (Bottier et al. 2017; Choudhury et al. 2023). Note that the magnitude of the $c_0$ coefficient is measured ≈ 90% of $W/(2P)$, the non-dimensional prefactor for velocity in the villi zone (see (A 8)). For the radial component, $d_0 \approx 0$ since the fluid is radially confined by walls and no net radial flow can occur.

We compare the modelled effective boundary velocity components ($\tilde{u}_{z,r}^e$) with the simulated instantaneous axial and radial velocities ($\tilde{u}_{z,r}$) at the villi tips ($r = 1.02H$), in figures 10 and 11, respectively. The comparison is made for four time-steps, for $\Delta\phi = \pi/2$ and $\pi/4$, for $\tilde{a} = 0.1$, using coefficient values identified from figure 12. The axial and radial velocity signals measured at $r \approx H$, which carry complicated perturbations at the scale of the villi, are described well by our simplified effective boundary velocity model. We also note that the model fit improves significantly for the smaller value of $\Delta\phi = \pi/4$. We remark here that the model in (5.3) and (5.4) implicitly assumes that $\Delta\phi = 2\pi/N$ can only decrease via an increase in the periodic domain size $L_z$ (i.e. an increase in the number of equispaced periodic villi $N$). The model will require modifications in the scenario where an increasing number of villi are packed within the same periodic distance $L_z$, by decreasing the intervillous space ($P-W$), and thus decreasing $\Delta\phi$.

The effective velocity model in (5.3) and (5.4) indicates a phase-lock between axial and radial velocity oscillations at the scale of $L_z$. This phase-lock is responsible for the unequal clockwise and counter-clockwise vortical flows within the mixing layer, as seen in figures 5 and 6. We confirm this by carrying out independent simulations with a planar (villi free) bottom boundary, where the axial and radial velocities are modelled only using the second terms in (5.3) and (5.4) (terms associated with $c_1$ and $d_1$). Furthermore, one can transform (5.3) and (5.4) in the reference frame of the propagating wave, by substituting $z^R = z + ct$. Under such transformation we immediately see that the perturbations occurring at the villi scale (terms associated with $c_2$ and $d_2$) will travel in the $+z$ direction, at the wave speed $c$. Our proposed model therefore captures the essential physics of villi scale non-reciprocal motion responsible for irreversible flow in the $+z$ direction in the Stokes regime.

To summarize, the effective boundary velocity model consists of three terms for each component. The first of these is a constant term responsible for irreversible flow, while the remaining two are propagating harmonic terms. The constant term is set to a non-zero value for the effective axial boundary velocity $\tilde{u}_z^e$ generating a net axial steady streaming flow in the channel while that for the radial $\tilde{u}_r^e$ effective boundary velocity is set to zero. For both components, the first of the harmonic terms imposes travelling contractions-relaxations at the longer wavelength of the periodic domain $L_z$. The second harmonic term adds villi scale perturbations, caused by the discrete nature of the villi-wall, to the long-wavelength boundary velocity wave.

## 6. Mixing boundary layer and fluid inertia

The classical Stokes second problem of an oscillating flat plate, bounding a semi-infinite fluid above it, predicts that a viscous boundary layer would manifest normal to the plate, named as the Stokes layer, scaling as $\delta_{St} \sim \sqrt{\mu/(\omega\rho)}$ (Schlichting 1960). In the Stokes flow regime, the fluid velocity simply follows that of the flat plate, and $\delta_{St} \to \infty$ ($Wo \to 0$). Stokes layer length scale $\delta_{St}$ decreases with faster oscillations (decreasing oscillatory to
<figure><img src="image_17.png" /><figcaption>Figure 13: Plot of the mixing layer height $l$ against the Womersley number $Wo = W/\delta_{St}$, for various $\Delta\phi$, for (A) $\bar{a} = 0.2$ and (B) $\bar{a} = 0.1$. When $Wo \ll 1$, we have viscosity dominant regime, and when $Wo \sim O(1)$ or greater we reach the fluid inertia dominant regime (shaded grey region). Note the plateauing of the curves when $Wo \ll 1$, and their rapid decrease with increasing $Wo$, when $Wo \gg 1$, indicating a flow regime transition. The measured mixing layer heights $l$ from (A) and (B), collapses as a double-power law when appropriately rescaled, as plotted in (C), delineating the two regimes. Note that in (C) we have not show the data-points for $\Delta\phi = 0$ and $\pi$. In (C), large x-axis shows the viscosity dominant regime, while at small x-axis shows the inertial regime. Contrast the scaling for $l$ in the inertial regime $(\propto \delta_{St}^{1/2})$ with that for the theoretical boundary layer $(\propto \delta_{St})$ for an oscillating flat plate (Schlichting 1960). Transition between the two regimes appears to be smooth, and we visually identify (0.83,6.66) as the critical transition point in (C). From this critical point, we obtain the functional transition curve: $l^c = 6.66e^{-\bar{a}W/Wo}$, which is plotted as dash-dotted (grey) lines in (A) and (B). These lines indicate the viscous to inertial regime transition, for all $\Delta\phi$ and $\bar{a}$.</figcaption></figure>

viscous time scale), and serves as a proxy for an inertial length scale, mapping the oscillatory fluid inertia with $Wo = W/\delta_{St}$. In our simulations we observe an emergent boundary layer in the instantaneous flow-field, called the mixing boundary layer (MBL), showing dissimilar dynamics to the classical Stokes layer. Curiously, the MBL height $l$ remains finite even in the Stokes flow regime (see figure 5), while the Stokes layer in this limit is large $\delta_{St} > R$. As discussed in section 5.1, the MBL is characterised by vortical flows, depends on $\Delta\phi$ and is further affected by the dynamic confinement of the Stokes layer.
## 6.1. Mixing boundary layer height

The MBL height $\ell$ separates the vortical flow region from a unidirectional flow region (advecting layer) above it. In our simulations, the axial velocity in the advecting layer is always positive, i.e. $u_z > 0$. We use this fact to characterize the height of MBL, by first decomposing the instantaneous axial velocity according to its direction as,

$$ u_z^{-}(r, z, t) = \begin{cases} u_z, & \text{if } u_z < 0 \\ 0, & \text{otherwise,} \end{cases} \quad \text{and} \quad u_z^{+}(r, z, t) = u_z - u_z^{-} \quad (6.1) $$

and then computing the time and axial average of the negative velocity $u_z^-$ as,

$$ \langle u_z^{ss-} \rangle_z(r) = \frac{1}{L_z T} \int_0^{L_z} \int_t^{t+T} u_z^{-}(t, z, r) dt dz, \quad (6.2) $$

which depends only on the radial distance $r$. We now identify the MBL height $\ell$ as the shortest radial distance measured from the villi base at which the $\langle u_z^{ss-} \rangle_z(r := \ell) = 0$. With this definition we obtain the MBL height $\ell$, measured from the villi base, within which the instantaneous velocity field is characterized by significant negative axial flow $u_z^-(z, r, t)$, and above which the $u_z^-(z, r, t)$ is negligible.

We measure $\ell$ for all simulated cases, and its non-dimensional value is plotted against increasing $Wo$ in figure 13 for various $\Delta\phi$, for (A) $\tilde{a} = 0.2$ and (B) $\tilde{a} = 0.1$. We first look at the non-propagating case with $\Delta\phi = 0$ (star marker) and $\Delta\phi = \pi$ (circle marker). For synchronised oscillations of the villi, $\Delta\phi = 0$, $\ell$ remains constant at the maximum radial distance of $r = R$ for small $Wo$, and then sharply decreases with a slope $\sim \delta_{St}$ at increased $Wo$. This linear decrease of $\ell$ with increased $Wo$, mimics the dynamics due to an oscillating flat plate, except at $Wo = 5.0$, where structure of the villi-wall boundary causes deviations. Here, the mixing boundary layer is controlled by the Stokes layer, and the effect of the villi-wall can be approximated by that of an oscillating flat plate positioned at an offset to the villi-wall base. Puthumana Melepattu & de Loubens (2022) have shown that the mixing layer for synchronous oscillations is entirely an inertial phenomenon and can be controlled by changing the intervillous confinement. On the other hand, when $\Delta\phi = \pi$, we have neighbour villi that oscillate exactly out of phase and generate a mixing layer. This mixing layer decays till the channel center-line $r = R$, and there is no advecting layer generated.

A distinct class of behaviour is seen in the MBL ($\ell$) plots for propagating contractions, $0 < \Delta\phi < \pi$. Here, $\ell$ plateaus to a constant value at low $Wo$, and is a function solely of $\Delta\phi$ and $\tilde{a}$, consistent with (2.5). The plateau value of $\ell$ decreases as $\Delta\phi$ increases from $\pi/5$ to $2\pi/3$. At high $Wo$ in the inertial regime, $\ell$ decreases with a slope $\sim \sqrt{\delta_{St}}$, contrasting with the linear scaling seen for $\Delta\phi = 0$. The slower relative decay of the MBL compared to $\delta_{St}$ indicates additional radial fluid pumping mechanisms due to intervillous contractions at play. Comparing figures 13(A) and (B), we also see that decreasing $\tilde{a}$ causes an increase in the relative levels of the measured MBL height. At smaller $\Delta\phi$, e.g. for $\Delta\phi = \pi/5$ (square marker), the Stokes to inertial transition of $\ell$ occurs at a smaller $Wo$, compared to e.g. $\Delta\phi = \pi/2$ (triangle marker). The critical transition $Wo$ is smaller for larger plateau values of $\ell$ (i.e. smaller $\Delta\phi$ values). Here, for $\Delta\phi \le \pi/4$, the increasing mixing layer height is truncated by the radial confinement of the geometry. For these cases, the mixing layer extends across the entirety of the lumen at $Wo \ll 1$, and the advecting layer disappears.

The plots also appear to show an apparent increase of $\ell$ at high $Wo = 5.0$, for certain cases, namely for ($\tilde{a} = 0.2, \Delta\phi = \pi/5$) (see supplementary movie 6), and for ($\tilde{a} = 0.1, \Delta\phi = \pi/5$) and ($\tilde{a} = 0.1, \Delta\phi = \pi/4$) (see supplementary movie 7). Upon further investigation we find that the flow structure in these cases is unlike that observed in the rest of the simulations. Here, vortical flow pattern emanating from the villi is altered, with the smaller counter-clockwise
<figure><img src="image_18.png" /><figcaption>Figure 14: Snapshots of the purely oscillatory velocity field, visualized by subtracting a constant velocity ($U_z^{AL}$, 0) from the instantaneous velocity $\mathbf{u}$, shown for two time-fractions $t/T = 0.45$ and $0.75$, for the case $\Delta\phi = \pi/2$, $Wo = 0.16$ and $\tilde{a} = 0.2$. $U_z^{AL}$, computed from (6.3), is a measure of the irreversible axial flow velocity in the advecting layer. The dashed (magenta) separation line between the mixing and advecting layers is reproduced from figure 5(A) for reference.</figcaption></figure>

vortical flow assuming a distinct hourglass-like shape, while the larger clockwise vortical flow reaching the channel center. Despite the high oscillatory inertia at $Wo = 5.0$, we have a mixing layer height that reaches the channel center. Here, the steady streaming flow gets significantly weakened (see figure 8, panel: $\Delta\phi = \pi/5$, $Wo = 5.0$) and the oscillatory flow dominates throughout the luminal space. This changed flow structure indicates the existence of another inertial flow transition beyond $Wo \gtrsim 5$, which would also depend on $\Delta\phi$ and $\tilde{a}$. Here the computed averaged MBL height ($\ell$) is inconsistent due to the change in the flow structure for these cases. In this paper we limit ourselves to $Wo \le 5.0$, and a well defined $\ell$, and therefore exclude these three points from further analysis in the present section.

The mixing layer height data shows a phenomenological collapse onto a double power-law master curve, as shown in figure 13(C). Since the data for non-propagating cases ($\Delta\phi = 0$ and $\pi$) exhibits a different class of behaviour, it does not collapse onto the master curve, and is therefore not shown. The $x$ variable for the collapse is $\delta_{St}\Delta\phi^{2/3}/H$ and the $y$ variable is $\ell/(e^{-\tilde{a}}\delta_{St})$. For large values of the $x$ variable, the flow is dominated by viscous effects, leading to the power law relation $y \propto 1/x$ (blue line), and thus $\ell \propto \Delta\phi^{-2/3}e^{-\tilde{a}}$. Here the mixing layer height $\ell$ is unaffected by oscillatory fluid inertia and therefore independent of $\delta_{St}$. On the other hand, as the $x$ variable decreases, the flow transitions into the oscillatory inertial regime, with the power law relation $y \propto x^{-1/2}$ (red line). The mixing layer height now depends on $Wo$ and scales as $\ell \propto \delta_{St}^{1/2}\Delta\phi^{-1/3}e^{-\tilde{a}}$. The double power law collapse therefore shows a regime transition from the Stokes regime (large $x$) into the inertial regime (small $x$), with a relatively smooth transition in between. At smallest values of $x$ variable, we see that the data collapse starts to fray, suggesting that this flow structure will not sustain at highly elevated oscillatory inertia. We visually identify the critical transition point between the two regimes as ($x \approx 0.83$, $y \approx 6.66$). This allows us to deduce a critical transition line in terms of $\Delta\phi$ as: $\ell^c = 5.5e^{-\tilde{a}}H\Delta\phi^{-2/3}$. And similarly in terms of $Wo$ as: $\ell^c = 6.66e^{-\tilde{a}}W/Wo$. This second critical transition line is overlaid over the unscaled $\ell/H$ data in figures 13(A) and (B) (grey, dash-dotted line), highlighting the onset of the oscillatory fluid inertia in the system. The inertial regime is shown by the grey shaded region in all plots in the figure 13.

We now look into the physical mechanism that leads to the appearance of distinct mixing and advecting layers in certain simulated cases, while not in others. Towards this end, we subtract the characteristic flow velocity $U_z^{AL}$ of the advecting layer from the entire flow-field.
<figure><img src="image_19.png" /><figcaption>Figure 15: Axially averaged $\langle\langle\cdot\rangle\rangle_z$ z-components of two important velocity measures, (A) the oscillatory axial velocity $\langle u_z^{os} \rangle_z$ (dashed lines) and (B) the steady streaming axial velocity $\langle u_z^{ss} \rangle_z$ (solid lines), plotted on the x-axis. These measures are computed as functions of radial distance $r/H$ (on y-axis), from equations (6.4) and (6.5) for various $\Delta\phi$ and plotted (column-wise) for three increasing Womersley numbers $Wo = 0.16, 1.58$ and $5.01$. The plots demonstrate that the oscillatory flow decays exponentially in $r$, while the steady flow (which is mostly uniform for low $Wo$) does not.</figcaption></figure>

Figure 14 shows this shifted velocity field for $\Delta\phi = \pi/2$ and $Wo = 0.16$, with the height of the MBL reproduced (dashed magenta line) from figure 5(A). Here, the axial component of the characteristic velocity is computed as,

$$U_z^{AL} = \frac{1}{HL_z} \int_0^{L_z} \int_{5H}^{R} u_z dr dz, \qquad (6.3)$$

making sure that we integrate well above the MBL height, and the radial component is set to $U_r^{AL} = 0$. The shifted velocity field in figure 14 comprises entirely of counter-rotating semi vertical flow structures that now extend till the central symmetry boundary. This structure is similar to the one seen when the mixing layer extends till the center line (e.g. figure 5(B)). This pattern is produced irrespective of $\Delta\phi$ in the Stokes flow regime ($Wo = 0.16$), demonstrating that the total flow-field is simply a sum of the irreversible steady streaming and oscillatory components.

The distinct mixing and advecting layers seen in the simulations manifest due to the competition between the strength of the oscillatory flow, which generates the counter-rotating vertical structures, and the irreversible steady flow shown in figure 4, along the radial direction. In order to understand the MBL height $\ell$ behaviour, we map the radial dependence
of the irreversible and vortical flows as,

$$ \langle u_z^{ss} \rangle_z(r) = \frac{1}{L_z} \int_0^{L_z} u_z^{ss}(z, r) dz \qquad (6.4) $$

$$ \langle u_z^{os} \rangle_z(r) = \frac{1}{2TL_z} \int_0^{L_z} \int_t^{t+T} |u_z(t, z, r) - u_z^{ss}(z, r)| dtdz, \qquad (6.5) $$

where, $\langle u_z^{ss} \rangle_z$ and $\langle u_z^{os} \rangle_z$ are the axially averaged steady streaming and oscillatory axial velocities. Figure 15 compares the strength of these two velocity measures, for three increasing values of $Wo$ (column-wise panels) and all $\Delta\phi$. In figure 15(A), the oscillatory velocity $\langle u_z^{os} \rangle_z$ shows an exponential decay with increasing $r/H$, for various $\Delta\phi$. This decay becomes more pronounced as $Wo$ increases. In contrast, as seen in figure 15(B), the steady streaming $\langle u_z^{ss} \rangle_z$ remains nearly constant, starting from the villi tips to the centre of the lumen at lower $Wo$, while becoming increasingly bulged just above the villi tips at higher $Wo$.

In summary, for a given small $\tilde{a}$, the non-reciprocal motion of the villi-wall generates irreversible or steady streaming axial flow. This flow is nearly uniform, and depends only on $\Delta\phi$ when viscosity dominates over oscillatory inertia. The propagating contraction-expansion kinematics of the villi-wall also generate axially oscillating flow, seen in the form of counter-rotating vortical structures. The strength of this oscillatory flow decays exponentially along the radial direction, with a decay rate dependent on $\Delta\phi$. The height of the MBL is now determined by a competition between these two effects. Within the mixing layer the oscillatory axial flow predominates, while in the advecting layer its strength falls below that of the steady streaming axial flow. Hence, in viscosity dominated cases, $\ell$ depends only on $\Delta\phi$. As oscillatory inertia increases, a decreasing Stokes layer $\delta_{St}$ enforces a dynamic radial confinement on both the steady and oscillatory flows, causing a depletion of the irreversible flux near the channel center, as well as a faster radial decay of the oscillatory axial flux. This in turn causes the MBL height to decrease non-trivially in the inertial regime, as a function of both $Wo$ and $\Delta\phi$, unlike that for a simple oscillating flat plate.

## 7. Micro-mixing and transport

We put forth the idea of an artificial villi-wall driven flow pumping and micro-mixing, using propagating contractions, similar to the adaptation of artificial ciliary arrays for microfluidic mixing and transport applications (den Toonder et al. 2008; Shields et al. 2010; Gu et al. 2020). Unlike the artificial villi arrays, the protruding villi-like structures can be rigid, which could prove advantageous in confined microfluidic geometries or for processing highly viscous fluids. The contraction-expansion wave along the villi-wall simultaneously generates axial as well as radial fluxes. While axial pumping can be split into an oscillatory and a steady streaming part, there cannot be any irreversible pumping along the radial direction because of the confining channel walls. In this section, we aim to understand how these flow patterns can be modulated by adjusting the two control parameters, ($\Delta\phi$, $Wo$).

### 7.1. Luminal axial flux

Axial fluid transport can be quantified by oscillatory and steady streaming axial fluxes. Similar to the calculation of the axial fluxes in the villi zone (in (5.1) and (5.2)), we compute steady streaming flux in the lumen zone ($H < r \le R$) as,

$$ Q_z^{ss,l} = \frac{1}{L_z} \int_H^R \int_0^{L_z} u_z^{ss}(z, r) dz dr \qquad (7.1) $$
<figure><img src="image_20.png" /><figcaption>Figure 16: Plots for (A) oscillatory axial ($Q_z^{os,l}$) and (B) steady-streaming axial ($Q_z^{ss,l}$) fluxes measured within the lumen zone ($H < r \le R$), plotted against the phase lag $\Delta\phi$, for increasing $Wo$. Filled markers (coloured) are for $\tilde{a} = 0.2$ and empty markers (black) for $\tilde{a} = 0.1$. Unlike the fluxes in the villi zone, rescaled data does not follow simple scaling. Note that in (B) $Q_z^{ss,l}$, shows non-monotonicity with increasing $Wo$, with the maximum irreversible steady flux occurring when $Wo \approx 2.82$.</figcaption></figure>

and the corresponding oscillatory flux as,

$$Q_z^{os,l} = \frac{1}{2TL_z} \int_H^R \int_0^{L_z} \int_t^{t+T} |u_z(t, z, r) - u_z^{ss}(z, r)| dtdzdr. \quad (7.2)$$

Figures 16(A) and (B) show the axial oscillatory and steady streaming fluxes in the lumen, respectively. For $\Delta\phi = 0$, the axial oscillatory flux ($Q_z^{os,l}$ in figure 16(A)) shows large variation due to fluid inertia. For $\Delta\phi = 0$, $Q_z^{os,l}$ decreases with increasing $Wo$, from a maximum value of $\sim U_0H$ at $Wo = 0.16$ in the Stokes regime, to near zero at $Wo = 5.0$. This behaviour is a direct consequence of the increasing oscillatory inertia that progressively limits the axial flux nearer the villi tips. For $\Delta\phi > 0$, the non-dimensional axial oscillatory flux is clustered within a narrow range, 0.2–0.4. The flux shows gradual and monotonic decrease with increasing $\Delta\phi$, for a given $Wo$. Moreover, the $Wo$ at which the largest flux is observed now depends on the parameter ($\Delta\phi, \tilde{a}$). Here, unlike for the case of the villi zone in figure 9(A), the data does not collapse onto a simple trigonometric curve.

The irreversible steady streaming luminal flux is shown as a function of $\Delta\phi$ in figure 16(B). In the viscous-dominated flow regime at low $Wo \ll 1$, the steady-state axial flux data, $Q_z^{ss,l}$, collapses onto the functional form $\sim \sin(\Delta\phi)$, analogous to the behavior observed in the villi zone (figure 9(B)). Extrapolating the geometric scaling seen in the villi zone to the lumen zone as, $W(R-H)(\sin(\Delta\phi))/(2PH)$, we see this (dotted red curve) scaling is slightly larger than the data-points for $Wo \le 0.28$. This is consistent with the effective velocity boundary constant coefficient $c_0 \approx 0.9W/(2P)$, measured in section 5.3.1. Thus, in the Stokes flow regime, we can approximate the total steady streaming axial flux in the channel as $Q_z^{ss} = Q_z^{ss,v} + Q_z^{ss,l} \approx (WR \sin(\Delta\phi))/(2PH)$. However with an increase in oscillatory inertia, the flux data no longer shows a functional collapse. Similar to its oscillatory counterpart, the flux here shows a non-monotonic behaviour with increasing $Wo$, stemming from fluid inertia. Note that $Q_z^{ss,l}$ at $Wo = 5.0$ is lower than that at $Wo = 2.82$, for all $\Delta\phi$. Though increasing $Wo$ increases the axial velocity of the fluid in the lumen, it also decreases the radial height to carry this fluid flux due to dynamic radial confinement of the
<figure><img src="image_21.png" /><figcaption>Figure 17: The oscillatory radial flux $Q_r^{os}$ per intervillous gap evaluated from (7.4) is plotted (on x-axis) against radial distance $r/H$ (on y-axis) for $\bar{a} = 0.2$ and three different $\Delta\phi = \pi, \pi/3$ and $\pi/5$. For each $\Delta\phi$, curves for four increasing Womersley numbers $Wo = 0.16$ (solid), 0.28 (dotted), 0.5 (dashed) and 5.0 (dash-dotted) are shown. Note that radial flux curves at low $Wo$ (for $Wo = 0.16$ and 0.28) coincide exactly. Minor deviations from these coincident curves are seen with the onset of oscillatory fluid inertia, when $Wo = 0.5$, and large deviations are seen at increased oscillatory inertia for $Wo = 5.0$.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 18: Plots for (A) the peak value of radial oscillatory flux ($Q_r^{os,max}$), and (B) the steady-streaming ($Q_r^{ss,max}$) fluxes, plotted against the phase lag $\Delta\phi$, for increasing $Wo$. Filled markers (coloured) are for $\bar{a} = 0.2$ and the empty markers (black) for $\bar{a} = 0.1$. In both plots, the peak flux data-points collapse solely as functions of $\Delta\phi$ (for $0 < \Delta\phi < \pi$), when rescaled appropriately. The scaling functions are shown with dashed (green) and solid (red) lines in (A) and (B), respectively.</figcaption></figure>

flow near the villi tips (see figure 4). The flux therefore reaches a maximum at $Wo \approx 2.8$ and falls beyond $Wo > 2.82$.
## 7.2. Radial pumping

The radial pumping of fluid in the channel cannot lead to an irreversible radial flow due to the confining symmetric top and bottom channel walls. Instead, all of the radially pumped fluid, due to intervillous contractions/expansions recirculates, forming counter-rotating semi-vortical flow structures, as seen in figures 5 and 6. These flow patterns are crucial for enhancing micro-mixing within the channel, especially in regions nearer the villi. In the physiological context of the intestine, this would enhance subsequent absorption of gut contents by the intestinal walls (Lentle & de Loubens 2015).

Radial pumping is quantified using radial oscillatory and steady streaming fluxes. Due to radial flow confinement, mass conservation will ensure that a simple integration of the velocity field along the z-axis yields zero net flux. We therefore define the steady streaming radial flux per intervillous gap as,

$$Q_r^{ss}(r) = \frac{1}{2N} \int_0^{L_z} |u_r^{ss}(z, r)| dz, \qquad (7.3)$$

Since the radial steady streaming flux is not irreversible, the equation (7.3) captures the strength of the vortical structures observed in the SSF. The corresponding oscillatory radial flux per intervillous gap is computed as,

$$Q_r^{os}(r) = \frac{1}{2TN} \int_0^{L_z} \int_t^{t+T} |u_r(t, z, r) - u_r^{ss}(z, r)| dtdz. \qquad (7.4)$$

Here, $N = 2\pi/\Delta\phi$ is the number of intervillous gaps along axial periodic length $L_z = NP$. Unlike in the axial case, the radial fluxes per intervillous gap ($Q_r^{os/ss}$) retain their radial dependence and have equal magnitude in $\pm r$ directions.

Figure 17 shows the non-dimensional oscillatory radial flux $Q_r^{os}(r)$ on the x-axis, as a function of the radial distance along y-axis, for three values of $\Delta\phi = \pi, \pi/3$ and $\pi/5$, at $\tilde{a} = 0.2$. The flux curves are plotted for four increasing values of $Wo = 0.16$ (solid), 0.28 (dotted), 0.5 (dashed) and 5.0 (dash-dotted), in order to demonstrate in effect of increasing oscillatory inertia on their radial behaviour. The radial flux rises to a peak value just outside the villi tips and then falls to zero along $r$ for all cases. Unsurprisingly, we see that the peak radial flux is generated when $\Delta\phi = \pi$, when the intervillous spaces undergo maximum strain for a given villi displacement $a$. With a decrease in $\Delta\phi$, the peak flux value drops. With decreasing $\Delta\phi$, the radial distance over which the flux drops to zero also increases, with the fastest drop seen for $\Delta\phi = \pi$, and the slowest for $\Delta\phi = \pi/5$. This is consistent with the behaviour of the mixing layer height $\ell$ with decreasing $\Delta\phi$.

In the lumen zone, when $Wo \ll 1$, (i.e. $Wo = 0.16$ and 0.28), the flux curves are coincidental for all $\Delta\phi$. This overlap demonstrates that in the absence of fluid inertia, the $Q_z^{os}$ and its radial dependence is solely controlled by the geometric motion of the villi-wall, ($\Delta\phi$). As $Wo$ increases ($Wo \approx 0.5$), slight deviations away from the coincidental curves appear along $r$, indicating an onset of inertial effects. In the oscillatory inertia dominated regime ($Wo = 5.0$), $Q_z^{os}$ curves decay much faster to zero along $+r$ in the lumen. This again is a consequence of decreasing $\delta_{St}$, and the radial flux in the lumen is now a function of both $\Delta\phi$ and $Wo$. However, in the villi zone the radial flux profiles are governed solely by the phase lag $\Delta\phi$, and these always overlap irrespective of $Wo$. Thus, similar to the axial fluxes, the oscillatory radial flux in the villi zone remains unaffected by oscillatory fluid inertia due to greater geometric flow confinement.

The peak value of oscillatory radial flux along $r$ is identified as $Q_z^{os,max}(r^{max}) = \max(Q_z^{os})$ and that for the steady streaming radial flux as $Q_z^{ss,max} = \max(Q_z^{ss})$. Figure 18(A) plots the rescaled oscillatory radial flux as a function of $\Delta\phi$, for increasing $Wo$. Here, we see a
remarkable data collapse onto the functional form ~ $(2/\pi) \sin(\Delta\phi/2)$, for all $W_o$ and $\tilde{a}$.
This simple scaling form can be derived by solving for the maximum strain-rate experienced
by an individual intervillous gap in the villi-wall system, the details of which are given in
section A.2. Minor deviations are seen for $\Delta\phi = 0$, where the flux does not remain identically
zero for increasing $W_o$. When $\Delta\phi = 0$, the radial oscillatory flux is generated solely through
an inertial mechanism, and increases linearly with $W_o$ (figure not shown). The peak radial
oscillatory flux generated by this inertial mechanism ($\Delta\phi = 0$) is completely negligible when
compared to that generated due to intervillous contractions ($\pi/5 \le \Delta\phi \le \pi$), by at least
an order of magnitude, for all $W_o$. The peak oscillatory radial flux therefore remains largely
unaffected by fluid inertia ($W_o$) for $\Delta\phi > 0$. Figure 18(B) plots the peak values of the steady
streaming radial flux as a function of $\Delta\phi$, for different $W_o$ and $\tilde{a}$. When rescaled with $U_0H\tilde{a}$,
similar to the steady streaming axial flux, all the data points for the propagating wave cases
$(0 < \Delta\phi < \pi)$ collapse onto a master curve, which is a simple trigonometric function solely
of $\Delta\phi$. The peak flux values for the non-propagating wave cases ($\Delta\phi = 0$ and $\pi$) do not show
a collapse and are affected by changes in $W_o$ and $\tilde{a}$.

7.3. Oscillatory radial flux and fluid inertia

Even though the peak radial flux values are more or less independent of $W_o$, oscillatory
inertia has an influence on the radial flow field. As seen in figure 17, increasing $W_o$ causes
a faster decay of the radial flux $Q_r^{os}(r)$ along $+r$, and will therefore influence luminal
micro-mixing significantly. We quantify the sharpness of the flux decay by fitting a simple
exponential model, within an appropriate truncated luminal region as, $Q_r^{os}(r^{low} \le r \le r^{high}) = b_0e^{-r/\lambda_r}$. Here, since we are interested in the flux decay beyond the peak flux value,
we take the lower radial limit for the fit model as, $r^{low} \approx (r^{max} + 0.2H)$, where $r^{max}$ is the
distance at which we find the maximum oscillatory radial flux $Q_r^{os}(r := r^{max}) = Q_r^{os,max}$. The
upper radial limit for the fit is found when the radial flux falls to (numerically approximated)
zero, $Q_r^{os}(r := r^{high}) = 0$. In this analysis, we first carry out cubic spline fitting of the radial
flux curves. We then analyze the first and second derivatives of the cubic splines to define
the radial position of the numerical zero, and therefore set the radial distance $r^{high}$. Since
flow is radially confined, fluid inertia effects on the radial fluxes are not obvious. We expect
to see a signature of fluid inertia through an estimation of the fitting parameter $\lambda_r$, which we
term as the radial-decay length constant.

For all simulations, we obtain the radial-decay length constant $\lambda_r$, which is plotted against $\Delta\phi$ in figure 19(A), for various $W_o$ and $\tilde{a}$. Here, we shift the $\lambda_r/H$ data points measured for $\Delta\phi = 0$ to $\Delta\phi = 2\pi$, in order to make them visible on the log-log scale. Strikingly, extracted $\lambda_r$ falls onto two distinct levels, both of which diverge linearly as $\lambda_r/H \propto \Delta\phi^{-1}$, without any additional rescaling. For a given $W_o$, the value of $\lambda_r$ is seen to jump from the lower level onto the upper level with an increase in $\Delta\phi$. The emergence of two distinct levels here indicates a sharp transition in the radial fluid flow regime.

It merits mention that, just as the MBL height $\ell$ serves to indicate the viscous to inertial
transition, with an identifiable critical transition point (figure 13), for the axial flow field,
$\lambda_r$ similarly tracks the viscous to inertial transition in the radial flow field. To visualize the
regimes and their transition, we define the level number $K = \lambda_r \Delta\phi/H$, which assigns a
scalar value to the level selected by $\lambda_r$. We plot the level number $K$ as a colour function of
$(\Delta\phi, W_o)$, for $\tilde{a} = 0.1$ and 0.2 in figures 19(B) and (C), respectively. The two phase-maps
are highly similar, indicating that the magnitude of the small amplitude oscillation has a
negligible effect on the oscillatory radial flux profiles. On these phase-maps, red indicates
viscosity dominated regime with a high $K$ value, while blue indicate inertial regime with a
low $K$ value, with the sharp transition between.

The regime transition zone is affected by both $W_o$ and $\Delta\phi$, and we approximate this as
<figure><img src="image_23.png" /><figcaption>Figure 19: (A) The radial-decay length constant $\lambda_r$ obtained from analysing the radial oscillatory flux $Q_r^{os}(r)$, plotted against $\Delta\phi$ for all values of $Wo$. Filled markers (coloured) are for $\tilde{a} = 0.2$ and empty markers (black) for $\tilde{a} = 0.1$. We see that $\lambda_r \propto 1/\Delta\phi$ and assumes two distinct levels. The two levels are quantified with a level number $K = \lambda_r \Delta\phi/H$, which is mapped over $(\Delta\phi, Wo)$ space, for (B) $\tilde{a} = 0.1$ and (C) $\tilde{a} = 0.2$. A jump in $K$ indicates the onset of oscillatory inertia on radial flow. This transition is the equivalent in both (A) and (B), and is approximately marked by the dotted (grey) line. The dotted lines are identical in both maps, and thus the transition seems independent of $\tilde{a}$. The dash-dotted (magenta) line plots the functional transition curve, $Wo = 1.21\Delta\phi^{2/3}W/H$ for reference, obtained from the critical transition point via an analysis of the mixing boundary layer in figure 13.</figcaption></figure>

a double-power law function $Wo \propto \Delta\phi^q$, which is marked with a dotted (grey) line, and is identical in both the phase-maps. And for both phase-maps, the power-law exponent of this transition function is seen to change around $(\Delta\phi \approx 2\pi/5, Wo \approx 0.67)$. Since $\lambda_r$ is an inverse linear function of $\Delta\phi$, we postulate that this change of exponent comes about because of the limiting effect of the channel height $R/H = 5.6$. In both figures 19(B) and (C), we also mark the critical transition line (dash-dotted, magenta) obtained earlier in section 6.1 from the critical transition point in figure 13(C) as: $Wo \approx 1.21\Delta\phi^{2/3}W/H$. Here, we re-emphasize the distinction that the dash-dotted critical line is obtained in section 6, exclusively through
<figure><img src="image_24.png" /><figcaption>Figure 20: Phase-maps for axial to radial pumping ratio, for (A) $\tilde{a} = 0.1$ and (B) $\tilde{a} = 0.2$, over the simulated $(\Delta\phi, W_o)$ space. The colour map shows the ratio of the total irreversible steady streaming axial $(Q_z^{ss} = Q_z^{ss,v} + Q_z^{ss,l})$ to maximum oscillatory radial $(Q_r^{os,max})$ flux. Note the distinct maximal axial pumping regions in both plots. The two maps become equivalent when $Q_z^{ss}$ is rescaled with $\tilde{a}$ (see figure 24).</figcaption></figure>

an analysis of the axial flow-field (by defining $u_z^-(t, z, r)$), whereas the dotted double-power transition function is obtained from phase-maps that are constructed solely by analysing the radial flow field (via $|u_r(t, z, r)|$) within the lumen. We therefore deduce that for a given $\Delta\phi$, oscillatory fluid inertia impacts radial flow field at a higher $W_o$, compared to its effect on the axial flow field. This disparity arises due to flow confining boundary walls in the radial direction.

### 7.4. Micro-mixing and axial pumping regimes

Adopting the villi-wall boundary with propagating intervillous contractions towards a microfluidic application would require understanding which control parameters can be tuned to modulate micro-mixing or axial pumping. To advance this understanding, we compute the ratio $Q_z^{ss}/Q_r^{os,max}$, which serves as a proxy quantifying the relative effectiveness of axial pumping compared to micro-mixing. Here, axial fluxes in the villi ($Q_z^{ss,v}$) and the lumen ($Q_z^{ss,l}$) are combined to give the total irreversible axial steady streaming flux $Q_z^{ss} = Q_z^{ss,v} + Q_z^{ss,l}$. The ratio is plotted as scalar contour of $(\Delta\phi, W_o)$ in figure 20 for both $\tilde{a} = 0.1$ and $0.2$. In either phase-map we identify two dominant regimes, one where $Q_z^{ss}/Q_r^{os,max} > 1.0$, for which the axial pumping dominates, whereas for $Q_z^{ss}/Q_r^{os,max} < 1.0$ we expect better relative micro-mixing in the lumen. Clearly, an increase in the oscillatory amplitude $\tilde{a}$, causes a greater portion of the sampled parameter space $(\Delta\phi, W_o)$ to prefer the axial pumping over micro-mixing, and vice versa. In both the maps, the ratio $Q_z^{ss}/Q_r^{os,max}$ is non-monotonous with a distinct maximum around $\pi/3 \le \Delta\phi \le \pi/2$, and $W_o \approx 2.8$, which should give best axial pumping relative to micro-mixing. Generally speaking, axial pumping is seen to dominate at high $W_o > 1.0$ in the system. Micro-mixing would do best at $\Delta\phi = \pi$ and $0.1 \le W_o \le 0.5$, and for low $W_o < 0.5$, its dominance falls gradually with decreasing $\Delta\phi$. Finally, we note that the two phase-maps become equivalent when the flux ratio is rescaled by $\tilde{a}$, i.e. $Q_z^{ss}/(\tilde{a}Q_r^{os,max})$ (see figure 24), and the role of the villi displacement amplitude is to simply shift the phase-map in the control parameter space $(\Delta\phi, W_o)$.
<figure><img src="image_25.png" /><figcaption>Figure 21: Flow-field at two time instances in the axially periodic domain of length $L_z \approx 23.4$ mm ($N = 73, P = 320$ µm), simulated with propagating effective boundary velocities at the bottom edge and symmetry at the top. The boundary wave travels from right to left, as seen from the progression of the streamlines with time. The colour field in (A) shows the instantaneous velocity magnitude while that in (B) shows the instantaneous shear strain rate in the fluid.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>Figure 22: Flow-field at two time instances in the axially periodic domain of length $L_z \approx 7.7$ mm ($N = 24, P = 320$ µm), simulated with propagating effective boundary velocities at the bottom edge and symmetry at the top. The boundary wave travels from right to left, as seen from the progression of the streamlines with time. The colour field in (A) shows the instantaneous velocity magnitude while that in (B) shows the instantaneous shear strain rate in the fluid.</figcaption></figure>

## 8. Propagating contractions along the rat duodenum

In this section we investigate the flow-field generated in the duodenum of the rat undergoing propagating longitudinal contraction at near-physiological conditions. Using image analysis, Lentle *et al.* (2012) have obtained spatio-temporal mapping of propagating longitudinal contraction-relaxations waves on the outer walls of living rat duodenum. Their experiments estimate the wave speed for propagating contractions at $\approx 13.6$ mm/s, though here a direct measurement of the wavelength has not been recorded due to size limitations of the extracted duodenal segments. From our villi-wall motility model, we estimate the periodic length of the propagating wave as $L_z = c/f \approx 23.4$ mm. As the intervillous mean distance is $\approx 320$ µm (see table 1), the number of villi in a periodic domain along the z-direction is estimated as $N \approx 73$. And the wave speed is translated into the model as a constant intervillous phase lag of $\Delta\phi = 2\pi/73$. We thus simulate a duodenal section of periodic length $L_z = 23.36$
mm, undergoing propagating contraction-relaxations using the effective velocity boundary conditions (5.3) and (5.4), proposed in section 5.3.1. Using resolved villi (as done in all previous sections) for simulating a domain with 73 periodic villi would entail significantly higher computational costs.

In their study, Lentle et al. (2012) also map non-propagating contractions from the rat duodenal walls, measuring their averaged axial wavelength to $7.7 \pm 4$ mm. Multiple previous studies have carried out fluid simulations over similarly sized rat intestinal sections, with non-propagating wall contraction-relaxations, both with (Zhang et al. 2020; Wang & Brasseur 2017; Wang et al. 2010) and without villi (de Loubens et al. 2013). These simulations consider the fluid to be Newtonian with a viscosity of $\mu = 1.0$ mPa·s. To facilitate comparison with literature, we also simulate propagating contractions-relaxations over this shorter periodic domain of $L_z = 7.68$ mm. Here, with $P = 320$ µm, we get $\Delta\phi \approx 2\pi/24$ for the effective velocity model, and a wave speed of $c \approx 4.45$ mm/s. Taking an average strain rate of $26\%\text{s}^{-1}$ and an oscillation frequency of $\approx 34.8$ cycles per minute ($f \approx 0.58$ 1/s) (Lentle et al. 2012), we estimate a maximum strain displacement of $2a \approx 71.7$ µm over an intervillous distance of $P \approx 320$ µm. We therefore set $\tilde{a} = 0.18$ in both these simulations. The coefficients for the effective velocity model $c_{0-2}$ and $d_{0-2}$ are then calculated from the heuristic fitting relations derived in the Stokes flow limit, discussed in 5.3.1.

The rat gut is modelled as a 2D rectangular domain where the velocity at the bottom boundary is given by the equations (5.3) and (5.4), the top boundary enforces flow symmetry along $r$, while the right and left boundaries are periodic inlet-outlets. Note that unlike in all previous simulations the bottom boundary starts at $r = H = 500$ µm, and the simulated domain height is $(R-H)$, since the effective velocity model excludes the villi zone. Here, we reduce the lattice resolution to $P = 16$ lattice nodes (= 320 µm), as the simulations no longer explicitly model the moving villi. In order to make comparison viable, we adopt $\mu = 1.0$ mPa·s and impose an oscillation frequency of $f = 0.58$ 1/s, which leads to a physiological Womersley number of $Wo = 0.38$. The periodic flow-field convergence is reached here with the same criteria as discussed in section 3.2. We present these simulations with the caveat that biologically realistic flow-fields in the intestine are expected to be strongly influenced by the non-Newtonian nature of the fluid, which is mainly constituted by a mixture of the ingested contents and the viscoelastic intestinal mucus. Furthermore, experimental evidence suggests that longitudinal muscle contractions are not radially symmetric, and also show large amplitude variations, both axially and temporally. We aim for our effective velocity model to inspire further detailed investigations into these complex and biologically relevant phenomena.

The instantaneous flow-fields for the two simulated cases with $L_z \approx 23.4$ mm and $L_z \approx 7.7$ mm are presented in figures 21 and 22, respectively. We show all results in this section in their respective physical units keeping in mind their bio-physiological implications. The flow-field streamlines are plotted at two time instances in panels (A) at $t/T = 0.45$ and in (B) at $t/T = 0.75$. In both the figures, the colour field in panel (A) plots the magnitude of the velocity field $|u|$, while that in (B) plots the fluid shear strain rate $\dot{\gamma}$. The shear strain rate is measured locally from the LBM simulation as $\dot{\gamma} = \sigma_{zr}/(2\mu)$, where $\sigma_{zr}$ is the viscous shear stress computed at each lattice node.

In both figures, the streamlines show two counter-rotating vortices, of approximate axial size $L_z/2$ and height $R-H$. These pair vortices advance along $-z$ due to applied wave velocity boundary condition. We see that the vortical structures (the mixing layer) extends till the channel center, indicating that we are in the micro-mixing dominated regime. This is unsurprising, and a consequence of the assumption used for the estimation of boundary velocity coefficients $c_{0-2}$ and $d_{0-2}$. In the case of the rat duodenum, the ratio $c_1/c_0 \propto 1/(\tilde{a} \sin(\Delta\phi)) \approx 249$, from coefficient relations given in 5.3.1, for $\Delta\phi = 2\pi/73$. This means
<figure><img src="image_27.png" /><figcaption>Figure 23: Plots of velocity measures for the rat duodenal simulations at physiological $Wo = 0.38$. Dashed (blue) and solid (red) curves are for slower ($c \approx 4.45$) and faster ($c \approx 13.55$) pendular-waves over periodic domains $L_z \approx 7.7$ mm and $L_z \approx 23.4$ mm, respectively. $\langle u_z^{ss} \rangle_z$ is computed from (6.4), $\langle u_z^{os} \rangle_z$, and similarly the radial component $\langle u_r^{os} \rangle_z$ from (6.5), and $\langle |\dot{\gamma}| \rangle_{t,z}$ from (8.1).</figcaption></figure>

that the oscillatory axial flow dominates the irreversible SSF all throughout the gut. We expect this estimation of the model coefficients to be improved upon, with future experiments that can directly map the motility along the inner villi-lined walls of the rat duodenum.

Interestingly, we also see that the vortex pair in both cases appears skewed. This happens because the background SSF is affected by fluid inertia. Nearer the channel center $r = R$, we observe greater steady axial flow velocities $u_z^{ss}$ relative to those near the bottom boundary (see figure 23). This inertial effect even at the relatively low $Wo = 0.38$ gives the pair-vortices their skewed character. We also observe that the effective velocity boundary condition on the bottom wall imparts the effect of discrete villi on the flow-field, which is especially visible on the shear strain rate in panel (B). This effect is seen to be more pronounced in the shorter domain ($L_z \approx 7.7$, figure 22) simulation, as the intervillous pitch becomes relatively more important. The maximum velocity magnitude in both simulated cases is nearly equal at ~100 µm/m. The maximum shear strain rate values and velocity magnitudes are seen to occur near the bottom boundary, as expected.

We quantify radial dependence of the velocity field and shear rate as follows. We compute the axially averaged steady streaming axial velocity $\langle u_z^{ss} \rangle_z$ using (6.4). We use (6.5) to compute the axially averaged oscillatory axial velocity $\langle u_z^{os} \rangle_z$, and similarly the corresponding oscillatory radial velocity $\langle u_r^{os} \rangle_z$. Further, we compute the axially averaged shear strain rate magnitude as,

$$ \langle |\dot{\gamma}| \rangle_{t,z}(r) = \frac{1}{TL_z} \int_{0}^{L_z} \int_{t}^{t+T} |\dot{\gamma}| dtdz. \qquad (8.1) $$

We plot the radial variations of these averaged flow-field measures in figure 23, with solid line for the case $\Delta\phi = 2\pi/73$ ($L_z \approx 23.4$ mm) and dashed line for $\Delta\phi = 2\pi/24$ ($L_z \approx 7.7$ mm). We see that the magnitude of the steady streaming velocity is at least an order of magnitude smaller than that of the oscillatory velocity components, in the respective simulations. The steady streaming velocity is seen to increase nearer the channel center than at the bottom wall due to inertial effects, which consequently causes a skewness of the vortex pair in the instantaneous flow. Examining, the oscillatory velocity measures we see that $\langle u_z^{os} \rangle_z$ has a minimum, and $\langle u_r^{os} \rangle_z$ has a maximum, at approximately 500 µm above the villi tips. We postulate that this is an indicator of the micro-mixing zone seen in Wang & Brasseur (2017), and would significantly enhance mixing and absorption nearer the villi. For $\Delta\phi = 2\pi/73$, the
maximum shear strain rate at the wall is $\dot{\gamma}^{\text{max}} \approx 0.21 \text{ 1/s}$, and the shear rate $\langle |\dot{\gamma}| \rangle_{t,z}$ shows an exponential decay from the villi tips till reaching the channel center line. However, the shorter domain at $\Delta\phi = 2\pi/24$, shows larger maximum wall shear strain rate of $\dot{\gamma}^{\text{max}} \approx 0.40 \text{ 1/s}$, but the shape of its decay to the channel center becomes more complex. Using a typical oscillatory axial velocity of 10 $\mu$m/s just above the villi tips, we estimate a Péclet number of $O(10)$ for a glucose molecule, while that for a typical gut bacterium such as *lactobacillus* would be of $O(10^3)$. Hence we expect the nature of the flow-field to play an important role in digestion and absorption of nutrients in the small intestines.

## 9. Conclusions

We investigate the fluid transport and mixing induced by propagating longitudinal contractions in an axially periodic channel — a prototypical motility pattern termed as ‘pendular-wave motion’, observed along the walls of the small intestine in various species, including humans. Beyond its physiological relevance, this mechanism also has potential applications in microfluidics, as it enables both pumping and mixing through a travelling contraction-expansion wave along the patterned channel wall. The channel walls consists of rigid structures of height $H$, called villi, protruding into the channel lumen. These villi have non-zero width ($W$) and are separated by intervillous gaps with an average axial mean distance $P$. In this study, the geometry of the villi-patterned channel wall is inspired from the leaf-like (or ridge-like) villi of the rat duodenum.

We use the two-relaxation-time lattice Boltzmann method (TRT-LBM) with the incompressible equilibrium to numerically solve the fluid flow in the intervillous gaps and the lumen of a 2D mirror symmetric channel. For describing the moving villi-wall boundary, we adopt a computationally local interpolated bounce-back scheme, and an iterative fresh-node refill procedure. The approach enables stable simulations with accurate capture of second-order flow phenomena, while preserving the locality and fidelity of the method.

Each villus along the wall oscillates harmonically with circular frequency $\omega$, with a constant phase lag $\Delta\phi$ relative to its neighbours. This phase lag governs both the wave speed $c$ and the axial domain length $L_z$, ensuring axial periodicity of the villi-wall motion. The villi oscillation frequency is non-dimensionalized as the Womersley number $Wo = W/\sqrt{v/\omega}$, and the oscillation amplitude as $\tilde{a} = a/W$.

Our simulations reveal an emergent boundary layer in the instantaneous flow-field, separating the flow into two horizontal layers. The layer immediately above the villi-wall is characterised by counter-rotating, asymmetrical semi-vortical flow structures, and we term this as the *mixing layer*. Intervillous contractions pump fluid in and out into the channel lumen causing re-circulatory flow between nodes and anti-nodes of the travelling wave along the channel boundary, thereby manifesting the mixing boundary layer. Above this layer we observe nearly uniform unidirectional axial flow, always in a direction opposite to the travelling wave along the villi-wall, which we term as the *advecting layer*. Remarkably, the two distinct layers persist even in the Stokes flow regime, $Wo \ll 1$.

Furthermore, we map the time-averaged flow-field $u^{ss}$, called the steady streaming flow (SSF), for various $(\Delta\phi, Wo)$ and small oscillation amplitude $\tilde{a} = 0.1$ and $0.2$. We see that the propagating villi-wall motion produces irreversible axial SSF throughout the channel, in a direction opposite to that of the propagating wave. The SSF sustains even in the absence of fluid inertia, i.e. when $Wo \ll 1$. We demonstrate that the physical mechanism for the time-symmetry breaking in the Stokes flow regime is the same as that for a simple interlinked three-sphere 1D swimmer (Najafi 2004), and is enabled by non-zero finite $\Delta\phi$. Here, irreversible axial flow is generated due to coordinated motion of the villi (and the contracting intervillous spaces), tracing a non-reciprocal path, even though each individual villus moves reciprocally.
Using simple geometric arguments we provide a scaling for the irreversible (steady streaming) axial flux, which matches well with the flux measured from the simulations, especially within the highly confined villi zone. With an increase in oscillatory inertia, the steady flow get increasingly confined to a region adjacent to the villi tips.

We introduce an effective velocity boundary condition $\mathbf{u}^e(r \approx H)$ that enables the computation of the flow field at larger scales without explicitly modelling individual villi. This boundary condition is derived from a Fourier analysis of the flow velocity above the villi tips in the Stokes limit. Notably, this analysis reveals that the axial velocity component can be decomposed into three terms: a constant term, which drives the irreversible SSF, and two propagating harmonic terms. The first corresponds to a travelling contraction-expansion wave at the domain length scale $L_z$, while the second introduces travelling perturbations in the axial velocity at the scale of the intervillous distance $P$. A similar decomposition holds for the radial velocity component, except that its constant irreversible term vanishes due to mass conservation imposed by the confining walls in the radial direction.

We establish that the mixing boundary layer in the instantaneous flow-field emerges due to a competition between the steady streaming axial flow and the oscillatory axial flow, away from the villi tips. The axial flow in the mixing layer is dominated by the oscillatory component of axial velocity, which decays exponentially normal to the villi-wall. Alternatively, the advecting layer is dominated by the steady component of axial velocity, controlled by the coordinated motion of the villi. In the viscosity dominated regime (low $Wo$), the steady streaming axial flow is nearly uniform in the channel. Moreover, with increasing $Wo$, both the oscillatory and steady components of axial flow experience normal confinement towards the villi tips, thus affecting the mixing boundary layer height $\ell$ in a non-trivial manner. The increasing dynamic confinement is caused by a decreasing viscous penetration layer, classically quantified by a length scale called the Stokes layer $\delta_{St} = \sqrt{v/\omega}$.

We achieve a phenomenological collapse for the mixing boundary layer height $\ell$ as a double power-law function. The collapse helps us to delineate the viscosity dominated regime where $\ell \propto \Delta\phi^{-2/3}e^{-\tilde{a}}$, from the inertial flow regime where $\ell \propto \delta_{St}^{1/2} \Delta\phi^{-1/3}e^{-\tilde{a}}$. Hence, $\ell$ increases with decreasing $\Delta\phi$, and is ultimately limited by the channel radius $R$. The point where power-law exponent changes helps us identify a critical transition line for axial flow: $Wo \approx 1.21\Delta\phi^{2/3}W/H$, from viscous to inertial regime.

Radial pumping is quantified using absolute radial fluxes derived from both instantaneous and steady streaming (SSF) velocity fields. The oscillatory component of the radial flux attains its maximum amplitude just above the villi and subsequently exhibits an exponential decay along the channel radius. Remarkably, the peak flux values collapse onto simple trigonometric functions of the phase lag $\Delta\phi$. Thus, peak radial fluxes are related solely to the geometric motion of the villi-wall. We extract the radial-decay length scale $\lambda_r$ for the radial fluxes and build a phase-map visualizing a transition for the radial flow from viscous to inertial regime. In comparison to the axial flow critical transition line, we find that the viscous to inertial transition for the radial flow is delayed in the $(\Delta\phi, Wo)$ space due to flow confinement by the channel walls.

We compare the axial steady streaming (irreversible) to radial oscillatory (reversible) pumping by the villi-wall in the channel, to identify regions of either dominant micro-mixing or of axial advection. These phase-maps have a pronounced maximum, suggesting an optimal region in the parameter space, $(\pi/3 \le \Delta\phi \le \pi/2, Wo \approx 2.8)$, for maximizing steady axial flow pumping.

Finally, returning briefly to the rat duodenum, we demonstrate the utility of our effective velocity boundary model to estimate the flow patterns therein. We employ the villi-wall effective velocity boundary condition $\mathbf{u}^e|_{r \approx H}$ to solve for the flow at larger duodenal scales,
without explicitly resolving the flow around individual villi. Using available physiological data, we estimate that the flow in the rat duodenum operates in a micro-mixing dominated regime, suggesting that the primary function of these ‘pendular-wave’ contractions is to enhance *in-situ* mixing rather than augmenting propulsion. The flow-field exhibits counter-rotating vortical structures, with their axial extent directly proportional to the wavelength of the contraction-relaxation wave. The shape of these vortices is modulated by the effects of mild fluid inertia within the channel.

A more detailed experimental characterization of motility on the internal gut wall is needed to better understand its role in intestinal transport. This includes studies at both the scale of individual villi and the scale of the propagating wavelength. Additionally, the intestinal fluids, such as digesta and mucus, exhibit complex non-Newtonian properties, which could significantly alter the flow dynamics and requires further investigation.

Beyond its physiological relevance, this study suggests that bio-inspired ‘pendular-wave’ villi-wall motion could be leveraged in microfluidic applications. Similar to artificial cilia arrays, such a mechanism could enable efficient mixing, pumping, or flow control in confined microfluidic channels, opening avenues for innovative biomimetic designs.

**Supplementary data.** Supplementary movies are available.

**Acknowledgements.** RV thanks Irina Ginzburg for her advise on LBM advanced boundary condition, as well as Jérémy O’Byrne and Anjishnu Choudhury for helpful discussions on irreversibility in the Stokes flow regime.

**Funding.** Most of the computations presented in this paper were performed using the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is supported by Grenoble research communities. LRP is part of the LabEx Tec21 (ANR-11-LABX-0030) and of the PolyNat Carnot Institute (ANR-11-CARN-007-01). The authors thank Agence Nationale de la Recherche for its financial support of the project TransportGut, ANR-21-CE45-0015.

**Declaration of interests.** The authors report no conflict of interest.

**Author ORCIDs.** R. Vernekar, https://orcid.org/0000-0002-3166-6564; C. de Loubens, https://orcid.org/0000-0002-4988-9168; C. Loverdo, https://orcid.org/0000-0002-0888-1717; Martin Garic, https://orcid.org/0000-0002-5871-3754; Dácil Idaira Yánez Martín, https://orcid.org/0009-0005-5642-0033;

# Appendix A. Geometric scaling of radial and axial fluxes

Let the size of any intervillous gap between two adjacent villi along $\hat{z}$ be given as,

$$ \zeta_i = X_{i+1} - X_i \qquad (A 1) $$

Substituting from (2.2) and simplifying we have,

$$ \zeta_i = P - W + 2a \sin\left(\frac{\Delta\phi}{2}\right) \sin\left(\omega t + \frac{\Delta\phi}{2}(2i - 1)\right), \qquad (A 2) $$

and the expansion (or contraction) rate of the gap is,

$$ \dot{\zeta}_i = 2a\omega \sin\left(\frac{\Delta\phi}{2}\right) \cos\left(\omega t + \frac{\Delta\phi}{2}(2i - 1)\right). \qquad (A 3) $$

## A.1. Scaling irreversible axial flux in villi zone

The velocity at the center of the gap bounded by the $i^{th}$ and $(i + 1)^{th}$ villi using linear interpolation, and further trigonometric simplification is,

$$ V_i = \frac{U_i + U_{i+1}}{2} = \omega a \cos\left(\frac{\Delta\phi}{2}\right) \sin\left(\omega t + \frac{\Delta\phi}{2}(2i - 1)\right) \qquad (A 4) $$
<figure><img src="image_28.png" /><figcaption>Figure 24: Equivalence of phase-maps for axial to radial pumping ratio, when rescaled by $\tilde{a}$. The map shows contour lines (dashed) for cases with $\tilde{a} = 0.1$ and filled contours for cases with $\tilde{a} = 0.2$, plotted over the $(\Delta\phi, W_o)$ parameter space.</figcaption></figure>

Now, since the axial flux through each gap would be proportional to the gap size and fluid velocity within, we compute the velocity weighted by the gap size as,

$$ M_{z,i} = \zeta_i V_i = (P - W)\omega a \sin\left(\omega t + (2i - 1)\frac{\Delta\phi}{2}\right) \cos\left(\frac{\Delta\phi}{2}\right) \\ + a^2 \omega \sin^2\left(\omega t + (2i - 1)\frac{\Delta\phi}{2}\right) \sin(\Delta\phi). \quad (A 5) $$

Taking the time integral for the irreversible component, only the second term in (A 5) survives, and can be simplified to,

$$ \overline{M}_z = \frac{1}{T} \int_t^{t+T} M_{z,i} dt = \frac{a^2 \omega}{2} \sin(\Delta\phi) \quad (A 6) $$

The scaling for the steady streaming irreversible flux through any gap within the villi zone can be now be written as,

$$ Q_z^{ss,v} = \frac{HN}{L_z} \overline{M}_z = \frac{a^2 \omega H}{2P} \sin(\Delta\phi). \quad (A 7) $$

And when non-dimensionalized with characteristic villi flux ($U_0H = a\omega H$), this is expressed as,

$$ \tilde{Q}_z^{ss,v} = \frac{\tilde{a}W}{2P} \sin(\Delta\phi). \quad (A 8) $$

## A.2. Scaling peak oscillatory radial flux

The instantaneous flux driven by the contraction/expansion of any intervillous gap is $|\dot{\zeta}_i|H$. We postulate that the peak oscillatory radial flux per intervillous gap would scale as,

$$ Q_r^{os,max} = \frac{1}{2T} \int_t^{t+T} |\dot{\zeta}_i| H dt = \frac{a\omega H}{T} \sin\left(\frac{\Delta\phi}{2}\right) \int_t^{t+T} \left|\cos\left(\omega t + \frac{\Delta\phi}{2}(2i - 1)\right)\right| dt, \quad (A 9) $$
where, we make use of the fact that $sin(\Delta\phi/2) > 0$ for $0 \le \Delta\phi \le \pi$. Now, the integral can be simplified to,

$$Q_r^{os,max} = \frac{2a\omega H}{T} \sin\left(\frac{\Delta\phi}{2}\right) \int_{-T/4}^{T/4} |\cos(\omega t)| dt = \frac{2a\omega H}{\pi} \sin\left(\frac{\Delta\phi}{2}\right), \quad (A 10)$$

which gives the non-dimensional oscillatory radial flux as,

$$\tilde{Q}_r^{os,max} = \frac{2}{\pi} \sin\left(\frac{\Delta\phi}{2}\right), \quad (A 11)$$

## Appendix B. Relative axial pumping to micro-mixing map

Here we demonstrate that the two phase-maps, for the axial irreversible to radial pumping ratio, plotted in figure 20, are equivalent when rescaled with $\tilde{a}$. The contour levels in the two phase-maps are rescaled and overlayed in the figure 24, with dashed contour lines plotted for the cases with $\tilde{a} = 0.1$ and the filled contours plotted for the cases with $\tilde{a} = 0.2$.

## REFERENCES

* AMEDZROVI AGBESI, RICHARD J & CHEVALIER, NICOLAS R 2022 Flow and mixing induced by single, colinear, and colliding contractile waves in the intestine. *Physical Review Fluids* 7 (4), 043101.

* BOLURIAAN, SAID & MORRIS, PHILIP J 2003 Acoustic streaming: from rayleigh to today. *International Journal of aeroacoustics* 2 (3), 255–292.

* BOTTIER, MATHIEU, FERNÁNDEZ, MARTA PEÑA, PELLE, GABRIEL, ISABEY, DANIEL, LOUIS, BRUNO, GROTBERG, JAMES B. & FILOCHE, MARCEL 2017 A new index for characterizing micro-bead motion in a flow induced by ciliary beating: Part II, modeling. *PLOS Computational Biology* 13 (7), e1005552.

* BRENNEN, CHRISTOPHER & WINET, HOWARD 1977 Fluid mechanics of propulsion by cilia and flagella. *Annual Review of Fluid Mechanics* 9 (1), 339–398.

* BUTTON, BRIAN, CAI, LI-HENG, EHRE, CAMILLE, KESIMER, MEHMET, HILL, DAVID B., SHEEHAN, JOHN K., BOUCHER, RICHARD C. & RUBINSTEIN, MICHAEL 2012 A Periciliary Brush Promotes the Lung Health by Separating the Mucus Layer from Airway Epithelia. *Science* 337 (6097), 937–941.

* CASSELBRANT, ANNA & HELANDER, HERBERT F. 2022 Asymmetric mucosal structure, mesenteric versus antimesenteric, in mouse, rat, and human small intestines. *Physiological Reports* 10 (24), e15547.

* CHEN, LI, YU, YANG, LU, JIANHUA & HOU, GUOXIANG 2014 A comparative study of lattice Boltzmann methods using bounce-back schemes and immersed boundary ones for flow acoustic problems. *International Journal for Numerical Methods in Fluids* 74 (6), 439–467.

* CHOUDHURY, ANJISHNU, FILOCHE, MARCEL, RIBE, NEIL M., GRENIER, NICOLAS & DIETZE, GEORG F. 2023 On the role of viscoelasticity in mucociliary clearance: A hydrodynamic continuum approach. *Journal of Fluid Mechanics* 971, A33.

* COSTALONGA, M, BRUNET, PHILIPPE & PEERHOSSAINI, H 2015 Low frequency vibration induced streaming in a hele-shaw cell. *Physics of Fluids* 27 (1), 013101.

* CUI, SONGYUAN, BHOSALE, YASHRAJ & GAZZOLA, MATTIA 2024 Three-dimensional soft streaming. *Journal of Fluid Mechanics* 979, A7.

* DAUPTAIN, ANTOINE, FAVIER, JULIEN & BOTTARO, ALESSANDRO 2008 Hydrodynamics of ciliary propulsion. *Journal of Fluids and Structures* 24 (8), 1156–1165.

* DING, YANG, NAWROTH, JANNA C., MCFALL-NGAI, MARGARET J. & KANSO, EVA 2014 Mixing and transport by ciliary carpets: A numerical study. *Journal of Fluid Mechanics* 743, 124–140.

* ELGETI, JENS & GOMPPER, GERHARD 2013 Emergence of metachronal waves in cilia arrays. *Proceedings of the National Academy of Sciences* 110 (12), 4470–4475.

* ESSER, FALK, MASSELTER, TOM & SPECK, THOMAS 2019 Silent Pumpers: A Comparative Topical Overview of the Peristaltic Pumping Principle in Living Nature, Engineering, and Biomimetics. *Advanced Intelligent Systems* 1 (2), 1900009.

* FISHMAN, AARON, ROSSITER, JONATHAN M, LEONTINI, JS & HOMER, ME 2022 Mixing in arrays of villi-like actuators. *Physics of Fluids* 34 (9), 094112.
GINZBURG, IRINA 2025 The Lattice Boltzmann Method with Deformable Boundary for Colonic Flow Due to Segmental Circular Constrictions. *Fluids* **10** (2), 22.

GINZBURG, IRINA, SILVA, GONCALO, MARSON, FRANCESCO, CHOPARD, BASTIEN & LATT, JONAS 2023 Unified directional parabolic-accurate lattice Boltzmann boundary schemes for grid-rotated narrow gaps and curved walls in creeping and inertial fluid flows. *Physical Review E* **107** (2), 025303.

GINZBURG, IRINA, VERHAEGHE, FREDERIK & D'HUMIERES, DOMINIQUE 2008 Two-relaxation-time Lattice Boltzmann scheme: About parametrization, velocity, pressure and mixed boundary conditions. *Communications in computational physics* **3** (2), 427–478.

GU, HONGRI, BOEHLER, QUENTIN, CUI, HAOYANG, SECCHI, ELEONORA, SAVORANA, GIOVANNI, DE MARCO, CARMELA, GERVASONI, SIMONE, PEYRON, QUENTIN, HUANG, TIAN-YUN, PANE, SALVADOR, HIRT, ANN M., AHMED, DANIEL & NELSON, BRADLEY J. 2020 Magnetic cilia carpets with programmable metachronal waves. *Nature Communications* **11** (1), 2637.

HALL, JON & CLARKE, NIGEL 2020 The mechanics of cilium beating: Quantifying the relationship between metachronal wavelength and fluid flow rate. *Journal of Fluid Mechanics* **891**, A20.

HE, XIAOYI & LUO, LI-SHI 1997a Lattice Boltzmann Model for the Incompressible Navier-Stokes Equation. *Journal of Statistical Physics* **88** (3), 927–944.

HE, XIAOYI & LUO, LI-SHI 1997b Theory of the lattice Boltzmann method: From the Boltzmann equation to the lattice Boltzmann equation. *Physical Review E* **56** (6), 6811–6817.

HOSOYAMADA, YASUE & SAKAI, TATSUO 2005 Structural and mechanical architecture of the intestinal villi and crypts in the rat intestine: integrative reevaluation from ultrastructural analysis. *Anatomy and embryology* **210** (1), 1–12.

ISHIKAWA, TAKUJI 2024 Fluid Dynamics of Squirmers and Ciliated Microorganisms. *Annual Review of Fluid Mechanics* **56** (Volume 56, 2024), 119–145.

JAFFRIN, MY & SHAPIRO, AH 1971 Peristaltic pumping. *Annual review of fluid mechanics* **3** (1), 13–37.

KHADERI, S. N., DEN TOONDER, J. M. J. & ONCK, P. R. 2011 Microfluidic propulsion by the metachronal beating of magnetic artificial cilia: A numerical analysis. *Journal of Fluid Mechanics* **688**, 44–65.

KOTAS, CHARLOTTE W., YODA, MINAMI & ROGERS, PETER H. 2007 Visualization of steady streaming near oscillating spheroids. *Experiments in Fluids* **42** (1), 111–121.

KRÜGER, TIMM, KUSUMAATMAJA, HALIM, KUZMIN, ALEXANDR, SHARDT, OREST, SILVA, GONCALO & VIGGEN, ERLEND MAGNUS 2017 *The lattice Boltzmann method*, vol. 10. Springer.

LAMMERS, WIM J. E. P. 2005 Spatial and temporal coupling between slow waves and pendular contractions. *American Journal of Physiology-Gastrointestinal and Liver Physiology* **289** (5), G898–G903.

LASER, DANIEL J & SANTIAGO, JUAN G 2004 A review of micropumps. *Journal of micromechanics and microengineering* **14** (6), R35.

LAUGA, ERIC 2011 Life around the scallop theorem. *Soft Matter* **7** (7), 3060–3065.

LENTLE, RG, JANSSEN, PWM, DE LOUBENS, C, LIM, YF, HULLS, C & CHAMBERS. P 2013 Mucosal microfolds augment mixing at the wall of the distal ileum of the brushtail possum. *Neurogastroenterology & Motility* **25** (11), 881–e700.

LENTLE, RG & DE LOUBENS, C 2015 A review of mixing and propulsion of chyme in the small intestine: fresh insights from new methods. *Journal of Comparative Physiology B* **185** (4), 369–387.

LENTLE, RG, DE LOUBENS, C, HULLS, C, JANSSEN, PWM, GOLDING, MD & CHAMBERS, JP 2012 A comparison of the organization of longitudinal and circular contractions during pendular and segmental activity in the duodenum of the rat and guinea pig. *Neurogastroenterology & Motility* **24** (7), 686–e298.

LEVITT, M. D., STROCCHI, A. & LEVITT, D. G. 1992 Human jejunal unstirred layer: Evidence for extremely efficient luminal stirring. *American Journal of Physiology-Gastrointestinal and Liver Physiology* **262** (3), G593–G596.

LIM, YF, DE LOUBENS, CLÉMENT, LOVE, RJ, LENTLE, RG & JANSSEN, PWM 2015 Flow and mixing by small intestine villi. *Food & function* **6** (6), 1787–1795.

LIM, YUEN FEUNG, LENTLE, ROGER G, JANSSEN, PATRICK WM, WILLIAMS, MARTIN AK, DE LOUBENS, CLEMENT, MANSEL, BRADLEY W & CHAMBERS, PAUL 2014 Determination of villous rigidity in the distal ileum of the possum (trichosurus vulpecula). *PloS one* **9** (6), e100140.

LOISEAU, ETIENNE, GSELL, SIMON, NOMMICK, AUDE, JOMARD, CHARLINE, GRAS, DELPHINE, CHANEZ, PASCAL, D'ORTONA, UMBERTO, KODJABACHIAN, LAURENT, FAVIER, JULIEN & VIALLAT, ANNIE 2020 Active mucus-cilia hydrodynamic coupling drives self-organization of human bronchial epithelium. *Nature Physics* **16** (11), 1158–1164.

DE LOUBENS, CLÉMENT, LENTLE, ROGER G, HULLS, CORRIN, JANSSEN, PATRICK WM, LOVE, RICHARD J & CHAMBERS, J PAUL 2014 Characterisation of mixing in the proximal duodenum of the rat during
longitudinal contractions and comparison with a fluid mechanical model based on spatiotemporal motility data. *PloS one* 9 (4), e95000.

DE LOUBENS, CLÉMENT, LENTLE, ROGER G, LOVE, RICHARD J, HULLS, CORRIN & JANSSEN, PATRICK WM 2013 Fluid mechanical consequences of pendular activity, segmentation and pyloric outflow in the proximal duodenum of the rat and the guinea pig. *Journal of the Royal Society Interface* 10 (83), 20130027.

LOUDON, CATHERINE & TORDESILLAS, ANTOINETTE 1998 The Use of the Dimensionless Womersley Number to Characterize the Unsteady Nature of Internal Flow. *Journal of Theoretical Biology* 191 (1), 63–78.

MACAGNO, E. O. & CHRISTENSEN, J. 1980 Fluid Mechanics of the Duodenum. *Annual Review of Fluid Mechanics* 12 (Volume 12, 1980), 139–158.

MAILMAN, D., WOMACK, W. A., KVIETYS, P. R. & GRANGER, D. N. 1990 Villous motility and unstirred water layers in canine intestine. *American Journal of Physiology-Gastrointestinal and Liver Physiology* 258 (2), G238–G246.

MARMOTTANT, PHILIPPE 2024 Large vortices from soft vibrations. *Journal of Fluid Mechanics* 986, F1.

MARSHALL, WALLACE F & KINTNER, CHRISTOPHER 2008 Cilia orientation and the fluid mechanics of development. *Current opinion in cell biology* 20 (1), 48–52.

MARSON, FRANCESCO, THORIMBERT, YANN, CHOPARD, BASTIEN, GINZBURG, IRINA & LATT, JONAS 2021 Enhanced single-node lattice Boltzmann boundary condition for fluid flows. *Physical Review E* 103 (5), 053308.

MELVILLE, J., MACAGNO, E. & CHRISTENSEN, J. 1975 Longitudinal contractions in the duodenum: their fluid-mechanical function. *American Journal of Physiology-Legacy Content Publisher: American Physiological Society*.

NAJAFI, ALI 2004 Simple swimmer at low Reynolds number: Three linked spheres. *Physical Review E* 69 (6).

NAJAFI, ALI & GOLESTANIAN, RAMIN 2005 Propulsion at low Reynolds number. *Journal of Physics: Condensed Matter* 17 (14), S1203.

NAMA, NITESH 2024 Acoustic streaming: insights across reynolds numbers. *Journal of Fluid Mechanics* 982, F1.

PANDEY, ANUPAM, CHEN, ZIH-YIN, YUK, JISOO, SUN, YUMING, ROH, CHRIS, TAKAGI, DAISUKE, LEE, SUNGYON & JUNG, SUNGHWAN 2023 Optimal free-surface pumping by an undulating carpet. *Nature communications* 14 (1), 7735.

PURCELL, EDWARD M 1977 Life at low reynolds number. *American journal of physics* 45 (1), 3–11.

PUTHUMANA MELEPATTU, MIDHUN & DE LOUBENS, CLÉMENT 2022 Steady streaming flow induced by active biological microstructures; application to small intestine villi. *Physics of Fluids* 34 (6), 061905.

RILEY, N. 2001 Steady Streaming. *Annu. Rev. Fluid Mech.* 33 (1), 43–65.

SATIR, PETER & CHRISTENSEN, SØREN TVORUP 2007 Overview of structure and function of mammalian cilia. *Annu. Rev. Physiol.* 69, 377–400.

SCHLICHTING, HERMANN 1960 *Boundary Layer Theory*. New York, McGraw-Hill.

SHAH, ALOK S, BEN-SHAHAR, YEHUDA, MONINGER, THOMAS O, KLINE, JOEL N & WELSH, MICHAEL J 2009 Motile cilia of human airway epithelia are chemosensory. *Science* 325 (5944), 1131–1134.

SHAPIRO, ASCHER H, JAFFRIN, MICHEL YVES & WEINBERG, STEVEN LOUIS 1969 Peristaltic pumping with long wavelengths at low reynolds number. *Journal of fluid mechanics* 37 (4), 799–825.

SHIELDS, AR, FISER, BL, EVANS, BA, FALVO, MR, WASHBURN, S & SUPERFINE, R 2010 Biomimetic cilia arrays generate simultaneous pumping and mixing regimes. *Proceedings of the National Academy of Sciences* 107 (36), 15670–15675.

SINNOTT, MATTHEW D, CLEARY, PAUL W & HARRISON, SIMON M 2017 Peristaltic transport of a particulate suspension in the small intestine. *Applied Mathematical Modelling* 44, 143–159.

STROCCHI, ALESSANDRA & LEVITT, MICHAEL D 1993 Role of villous surface area in absorption science versus religion. *Digestive diseases and sciences* 38 (3), 385–387.

TABATA, OSAMU, HIRASAWA, HIROSHI, AOKI, SHINTARO, YOSHIDA, RYO & KOKUFUTA, ETSUO 2002 Ciliary motion actuator using self-oscillating gel. *Sensors and Actuators A: Physical* 95 (2-3), 234–238.

TAKAGI, D & BALMFORTH, NJ 2011 Peristaltic pumping of viscous fluid in an elastic tube. *Journal of Fluid mechanics* 672, 196–218.

TAO, SHI, HU, JUNJIE & GUO, ZHAOLI 2016 An investigation on momentum exchange methods and refilling algorithms for lattice Boltzmann simulation of particulate flows. *Computers & Fluids* 133, 1–14.

TATSUNO, MASAKAZU 1973 Circulatory streaming around an oscillating circular cylinder at low reynolds numbers. *Journal of the Physical Society of Japan* 35 (3), 915–920.
DEN TOONDER, JAAP, BOS, FEMKE, BROER, DICK, FILIPPINI, LAURA, GILLIES, MURRAY, DE GOEDE, JUDITH, MOL, TITIE, REIJME, MIREILLE, TALEN, WIM, WILDERBEEK, HANS & OTHERS 2008 Artificial cilia for active micro-fluidic mixing. *Lab on a Chip* **8** (4), 533–541.

WANG, YANXING & BRASSEUR, JAMES G 2017 Three-dimensional mechanisms of macro-to-micro-scale transport and absorption enhancement by gut villi motions. *Physical Review E* **95** (6), 062412.

WANG, YANXING, BRASSEUR, JAMES G, BANCO, GINO G, WEBB, ANDREW G, AILIANI, AMIT C & NEUBERGER, THOMAS 2010 A multiscale lattice boltzmann model of macro-to micro-scale transport, with applications to gut function. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences* **368** (1921), 2863–2880.

WESTERGAARD, H., HOLTERMULLER, K. H. & DIETSCHY, J. M. 1986 Measurement of resistance of barriers to solute transport in vivo in rat jejunum. *American Journal of Physiology-Gastrointestinal and Liver Physiology* **250** (6), G727–G735.

WU, JUNRU 2018 Acoustic streaming and its applications. *Fluids* **3** (4), 108.

ZHANG, JUNFENG 2011 Lattice Boltzmann method for microfluidics: Models and applications. *Microfluidics and Nanofluidics* **10** (1), 1–28.

ZHANG, YANAN, WU, PENG, JEANTET, ROMAIN, DUPONT, DIDIER, DELAPLACE, GUILLAUME, CHEN, XIAO DONG & XIAO, JIE 2020 How motility can enhance mass transfer and absorption in the duodenum: Taking the structure of the villi into account. *Chemical Engineering Science* **213**, 115406.