# 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. Journal of Fluid Mechanics, 2025, 1024, pp.A9. 10.1017/jfm.2025.10873. hal-05382981

HAL Id: hal-05382981

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

Submitted on 26 Nov 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¹†, Faisal Ahmad¹,³, Martin Garic ², Dácil Idaira Yánez Martín ¹,³, Claude Loverdo², Stéphane Tanguy³ and Clément de Loubens¹‡

¹Univ. Grenoble Alpes, CNRS, Grenoble INP, LRP, 38000 Grenoble, France

²Sorbonne Université, CNRS, IBPS, JJP, 75005 Paris, France

³Univ. Grenoble Alpes, CNRS, UMR 5525, VetAgro Sup, Grenoble INP, TIMC, 38000 Grenoble, France

Inspired by small intestine motility, we investigate the flow induced by a propagating pendular-wave along the walls of a channel lined with rigid, villi-like microstructures. The villi undergo harmonic axial oscillations with a phase lag relative to their neighbours, generating travelling patterns of intervillous contraction.

Using two-dimensional lattice Boltzmann simulations, we resolve the flow within the villi zone and the lumen, sampling small to moderate Womersley numbers. We uncover a mixing boundary layer (MBL) just above the villi, composed of semi-vortical structures that travel with the imposed wave. In the lumen, an axial steady flow emerges, surprisingly oriented opposite to the wave propagation direction, contrary to canonical peristaltic flows. We attribute this flow reversal to the non-reciprocal trajectories of fluid trapped between adjacent villi, and derive a geometric scaling law that captures its magnitude in the Stokes regime.

The MBL thickness is found to depend solely on the wave kinematics given by intervillous phase lag in the low-inertia limit. Above a critical threshold, oscillatory inertia induces dynamic confinement, limiting the radial extent of the MBL and leading to non-monotonic behaviour of the axial steady flux.

We further develop an effective boundary condition at the villus tips, incorporating both steady and oscillatory components across relevant spatial scales. This framework enables coarse-grained simulations of intestinal flows without resolving individual villi.

Our results shed light on the interplay between active microstructure, pendular-wave and finite inertia in biological flows, and suggests new avenues for flow control in biomimetic and microfluidic systems.

## 1. Introduction

Inspired by the 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 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).

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

<footnote>‡ Email address for correspondence: clement.de-loubens@univ-grenoble-alpes.fr</footnote>
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 (Ding et al. 2014; Hall & Clarke 2020). At the scale of an individual cilium, the beat is inherently non-reciprocal, effectively breaking the time reversibility of low Reynolds number flows. At 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 300–1000 µm in height (Garic et al. 2025). 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).

Assuming that villi passively follow the dynamics of mucosal deformations induced by contractions of smooth muscle (Lentle et al. 2013), one of the simplest motility patterns to consider is the forcing of villi by *pendular-wave* activity. This 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 of height *H* as rigid microstructures (Lim et al. 2014) undergoing harmonic oscillations. The underlying *pendular-wave* is manifested by an array of villi oscillating with a constant phase lag $\Delta\phi$ between adjacent villi (see figure 1). The wavelength $L_z$ and the wave speed *c* are therefore defined by the angular frequency $\omega$ and the number of villi per wavelength, given by $N = 2\pi/\Delta\phi$ in this array, to enforce periodic boundary conditions.

Whereas this framework resembles the metachronal waves observed in ciliary systems dominated by viscous effects (Ding et al. 2014; Hall & Clarke 2020), the villi are at least an order of magnitude larger than cilia, and thus inertial effects are expected to significantly influence the flow field induced by the oscillatory activity of a villus array (Wang & Brasseur
2017; Puthumana Melepattu & de Loubens 2022). In flows driven by oscillatory boundaries, fluid inertia introduces an additional hydrodynamic length scale, the Stokes layer thickness, defined as $\delta_{St} = \sqrt{\nu/\omega}$, where $\nu$ is the kinematic viscosity (Schlichting 1960). This length scale directly influences both the location and strength of the steady streaming flow (SSF), an effect with important implications for mixing in microfluidic systems (Costalonga et al. 2015; Fishman et al. 2022). Taking the villus width $W$ as the characteristic geometric length scale, the Womersley number defined as $Wo = W/\delta_{St}$, is of order unity under physiological conditions, confirming a departure from the Stokes flow limit (see Table 1). Although SSF is typically associated with inertial effects, it can also arise under 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, can be interpreted as a form of SSF in the Stokes regime, although this classification is not commonly adopted in the literature (Riley 2001; Hall & Clarke 2020).

Our aim is to decipher the physical mechanisms that govern the spatio-temporal organization of the flow when villi-like microstructures are forced by *pendular-wave* activity. In Section 2, we define the problem and perform the non-dimensionalization. The lattice Boltzmann method (LBM), used to compute the flow fields, is described in Section 3, along with a high-accuracy treatment of the moving boundaries required to capture secondary flow phenomena such as SSF. Section 4 presents the flow generated by non-propagating villus oscillations, a situation in which SSF manifests through inertial effects, and whose magnitude is several orders lower than that of the instantaneous flow. In Section 5, we show that a phase-lagged travelling-wave villus motion produces net directional transport and the emergence of a near-wall mixing boundary layer (MBL). Section 6 explores the structure and evolution of the MBL over a range of Womersley numbers from small to moderate values. In Section 7, we explain the physical origin of the axial steady flow opposing the pendular-wave and propose a mechanistic scaling law that captures its magnitude and direction. Section 8 introduces an effective boundary condition at the villus tips that reproduces the key flow features without resolving individual villi. Finally, conclusions and perspectives are presented in Section 9.

# 2. Problem statement and non-dimensionalization

## 2.1. *Problem statement*

To explore the physical mechanisms underlying the spatio-temporal organization of the flow, we performed a systematic parametric study using a symmetric 2D planar channel of diameter $2R$, whose walls are patterned with moving ridge-like villus structures of uniform height $H$ and width $W$, as illustrated in figure 1. Previous 3D simulations of stationary contractions by Lim et al. (2015), with different 3D villus shapes, show only minor qualitative differences in the resulting flow fields.

The fluid flow is solved both in the intervillous spaces and in the main channel lumen. The villi are modelled as rigid structures. To validate this assumption, we carry out a simple order of magnitude calculation considering a single villus as a slender cantilever beam of length $L \approx 500$ µm and square cross-section of width $W \approx 200$ µm. The villus beam, made up of healthy intestinal extracellular matrix, is given to possess a Young's modulus of $\approx 2.9$ kPa (Creff et al. 2021). At a large chyme viscosity of 100 mPa·s, a maximum wall strain rate of $\approx 0.47$s⁻¹ was computed for physiologically mapped pendular motion in the rat gut (de Loubens et al. 2013). Using the slender beam theory, purely viscous stresses here would result in a maximum villus tip deflection of < 0.5 µm. This is < 0.1% of the total villus length, validating the rigidity assumption for the villi. The rigid villi condition is further supported by two experimental studies (Lim et al. 2014; Lentle et al. 2013), which demonstrate that villus motion is driven by contractions of the underlying muscularis mucosa
<figure><img src="image_2.png" /><figcaption>Figure 1: (A) Leaf-like and ridge-like villi of the small intestine of rat (duodenum and ileum). (B) The planar 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. (D) The same seven villi viewed from the top (z, t axes), showing the contractions (in red) and expansions (in hatched purple) of the inter-villus gaps propagating in the −z direction.</figcaption></figure>

and does not result from bending deformation under physiological fluid shear stresses. Note
that in the brushtail possum, smooth muscle contractions induce wrinkling of the mucosa
(Lentle et al. 2013). We neglect this phenomenon in the present model, as our ongoing
experiments in the rat duodenum have not shown such wrinkling, suggesting that it may be
species-dependent.

The spatial coordinates along the axial and radial dimensions of the channel are given by
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 purpose of analysis the domain is parted into two zones; the region patterned with
villi where 0 ≤ r ≤ H is termed the villi zone, and the channel space where H < r ≤ 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 Lz = 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}, \quad (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}, \quad (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 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 finite number of 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 f P / \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$. The bottom walls of the intervillous gaps are assigned instantaneous velocities that are linearly interpolated between those imposed by their bounding villi.

The fluid is assumed Newtonian, with dynamic viscosity $\mu$ and mass density $\rho$. The intestinal fluids are generally highly heterogeneous, and the Newtonian assumption provides a well-defined baseline in the regime of low particulate volume fraction, particularly in the proximal part of the small intestine where the digesta is mixed with a large quantity of secretions (Hardacre et al. 2018). This assumption, which is also relevant in view of potential microfluidic applications, facilitates the elucidation of the underlying physical mechanisms and the comparison across studied parameters. The flow is modelled by the incompressible Navier-Stokes equations given by,

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

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

where $p$ is the pressure and $\mathbf{u}$ the fluid velocity vector field.

## 2.2. Non-dimensionalization

One observes that there are multiple length and time scales in the problem, as seen in table 1. We can identify time scales for four physical processes: (i) the villi oscillation timescale $t^\omega = 1/\omega$, (ii) the fluid advective timescale $t^{adv} = W/(\omega a)$, (iii) the viscous timescale $t^\mu = W^2\rho/\mu$, and (iv) the time scale imposed by the boundary travelling wave $t^c = W/c = W\Delta\phi/(\omega P)$.

Furthermore, oscillatory boundaries introduce a characteristic length scale for viscous dissipation normal to the wall, termed as the Stokes layer, $\delta_{St} = \sqrt{\mu/(\omega\rho)}$ (Schlichting 1960). A priori, it is unclear which length and time scales would dominate the problem at hand. 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). We make the choice for non-dimensionalization as: $\tilde{u} = \mathbf{u}/(\omega a)$, $\tilde{p} = pW/(\mu\omega a)$, $\tilde{x} = \mathbf{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 \tilde{\nabla} \tilde{u} = - \frac{1}{\tilde{a} W o^2} \tilde{\nabla} \tilde{p} + \frac{1}{\tilde{a} W o^2} \tilde{\nabla}^2 \tilde{u}. \quad (2.5)$$

Two non-dimensional numbers emerge in equation (2.5), namely the Womersley number $Wo = W/\delta_{St}$ (Loudon & Tordesillas 1998) and the reduced oscillation amplitude $\tilde{a} = a/W$.
<table><caption>Table 1: Typical values characterizing villus geometry in the small intestine (mouse, rat, possum, rabbit, chicken, human, horse; see details in Garic et al. 2025) and longitudinal motility in rat (Lentle et al. 2012; de Loubens et al. 2013), possum (Lentle et al. 2013), and rabbit (Lammers 2005; Fullard et al. 2014). Orders of magnitude of dimensionless parameters are estimated based on these physiological values.</caption><tr><td>Radius of the small intestine</td><td>R</td><td>3–30 mm</td><td>(Lammers 2005; Lentle et al. 2012, 2013)</td></tr><tr><td>Width of the villi</td><td>W</td><td>50 – 400 µm</td><td>(Garic et al. 2025)</td></tr><tr><td>Length of the villi</td><td>H</td><td>300 – 1000 µm</td><td>(Garic et al. 2025)</td></tr><tr><td>Pitch</td><td>P</td><td>30 – 480 µm</td><td>(Garic et al. 2025)</td></tr><tr><td>Frequency</td><td>f</td><td>0.1 – 0.7 Hz</td><td>(Lammers 2005; Lentle et al. 2012, 2013)</td></tr><tr><td>Wave speed</td><td>c</td><td>10 – 30 mm/s</td><td>(Lammers 2005; Lentle et al. 2012)</td></tr><tr><td>Velocity amplitude</td><td>U<sub>0</sub></td><td>0.2 – 5 mm/s</td><td>(de Loubens et al. 2013; Fullard et al. 2014)</td></tr><tr><td>Displacement amplitude</td><td>a = U<sub>0</sub>/ω</td><td>50 – 1000 µm</td><td>(de Loubens et al. 2013; Fullard et al. 2014)</td></tr><tr><td>Fluid density</td><td>ρ</td><td>~994 kg/m<sup>3</sup></td><td>water at 37°C</td></tr><tr><td>Fluid viscosity</td><td>μ</td><td>≥ 0.7 mPa.s</td><td>water at 37°C</td></tr><tr><td>Stokes layer</td><td>δ<sub>St</sub> = √μ/ωρ</td><td>400 – 1000 µm</td><td>estimated</td></tr><tr><td>Womersley number</td><td>Wo = W/δ<sub>St</sub></td><td>0.05 – 1.0</td><td>estimated</td></tr><tr><td>Phase lag</td><td>Δφ = ωP/c</td><td>10<sup>-5</sup> – 10<sup>-1</sup></td><td>estimated</td></tr></table>

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}. \qquad (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$.

The geometry of villi varies both within a single species and across different species, as listed in Table 1. However, the differences observed within the same region are of the same order of magnitude. Therefore, it is not relevant to study large variations in the geometric ratios. We thus fix these ratios to the typical values reported for the rat small intestine, i.e. $R/H = 5.6$, $H/W = 2.5$, and $P/W = 1.6$, for all simulations (Lentle et al. 2012; Hosoyamada & Sakai 2005; Casselbrant & Helander 2022).

The longitudinal muscular activity is reported at ~0.5 Hz in various species (table 1). Based on this value, the typical physiological Stokes layer thickness $\delta_{St}$ is estimated at least 470 µm for a watery digesta. Consequently, for typical rat villi lengths of ≈ 500 µm, the typical physiological $Wo = W/\delta_{St} \approx 0.42$. There is no direct measurement for $\Delta\phi$ reported in the literature; we estimate it to be small, with an upper bound of 0.2. In this paper we go beyond the limited physiological range of values of $Wo$ and $\Delta\phi$, in order to investigate the physical mechanisms governing the flow patterns and to assess their potential for microfluidic applications. $Wo$ is varied over three orders of magnitude and $\Delta\phi$ from 0 to $\pi$ in discrete steps. Regarding the amplitude of oscillations, we restrict our study to the small-amplitude regime, with $\tilde{a} = 0.1$ and $0.2$. This choice is motivated by previous studies on steady streaming flows, which have shown that more complex flow structures can emerge in the large-amplitude regime, even at low $Wo$ (Tatsuno 1973). When $Wo \ll 1$, along with $\tilde{a} \ll 1$, we see that the equation (2.5) reduces to the steady Stokes equation. However, when $Wo \gg 1$ we expect an inertial flow regime, with all terms in (2.5) playing a role.

The results presented in this paper are based on a total of 128 simulations spanning the parameter space ($Wo, \Delta\phi, \tilde{a}$).
### 3. Numerical methods

#### 3.1. Two-relaxation-time lattice Boltzmann solver

We solve the incompressible Navier-Stokes equations using the lattice Boltzmann method (LBM) using the D2Q9 lattice, which discretizes coordinate space in two directions and the velocity space in nine directions (He & Luo 1997b; Zhang 2011). The simulation code is validated for the Stokes second problem constrained to a symmetric channel to be second-order accurate with the comparisons shown in figure S1, and is open sourced (Vernekar et al. 2025). One key advantage of LBM is that it solves the flow on a fixed Cartesian grid, which remains efficient even in the presence of moving boundaries. This grid-based formulation allows for straightforward parallelization, as most operations are local. Moreover, 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_{\bar{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) units, 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 $g^{eq}$ from (3.4) for the next time step.

## 3.2. Moving boundary condition

Contrary to the LB simulations involving cilia (Hall & Clarke 2020), where the elongated microstructures are thinner than one lattice unit, villi have a finite thickness that spans several lattice nodes. 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 we take 20 lattice nodes for the width $W$ of each villi and 32 lattice nodes for the pitch $P$.

## 3.3. Fresh node treatment

The process of re-initializing 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^*, u^*)$ from (3.5).

(v) Then compute $g^{eq,*}(\rho^*, 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^*, u^*)$ (Marson et al. 2021). During the inner-iterations, we update $\rho, u$ and populations $g^{eq}, g$ only at the fresh nodes (and not at their link-wise neighbours).

### 3.4. Numerical convergence

The computations are initiated with the fluid at rest in the channel. As the simulation advances in time, the continuous evolution of the instantaneous flow field is obtained. At the end of every oscillation period $T$, we compute the developing SSF, $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,

$$u^{ss,*}(\mathbf{x})|_{T'} = \frac{\Delta t}{T} \sum_{i}^{T/\Delta t} 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 $L_2$ norm as,

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

When $E_{L2}$ falls below 1%, the SSF is considered to have converged, and the simulation is set to terminate after running for one additional time period, $t = T' + T$. During this final period, the time-periodic instantaneous velocity field $u(\mathbf{x}, t)$ and the steady streaming flow $u^{ss}(\mathbf{x})$ are recorded.

## 4. Non-propagating contractions and steady streaming flow

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 $\bar{a} = 0.2$, unless otherwise stated. Flow velocities are rescaled by the characteristic villus velocity $U_0 = \omega a$, and plotted streamline width is proportional to the local velocity magnitude to illustrate flow strength.

### 4.1. Instantaneous flow

Figure 2(A) shows the instantaneous flow-field for $\Delta\phi = 0$, for $Wo = 0.5$ at a given time instance. This is the 2D equivalent of the synchronised oscillations of an infinite array of villi recently studied in details by Puthumana Melepattu & de Loubens (2022). At low $Wo$, the flow closely follows the oscillations of the villus walls. As $Wo$ increases and $\delta_{St}$ decreases, inertial effects progressively localize the flow near the villus tips and introduce a lag in the lumen (see figure S2). This results in transient recirculation zones above the villus tips.

When $\Delta\phi = \pi$, neighbouring villi oscillate exactly out of phase, as shown in figure 2(B) for $Wo = 0.5$ at a given time step. This motion generates counter-rotating vortical flow in the domain. With increased $Wo$, these vortical structures localize near villi tips (similarly to $\Delta\phi = 0$ cases) due to radial confinement by a decreasing $\delta_{St}$ (see figure S3).
10

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

Figure 2: Flow induced by non-propagation oscillations of the villi over two neighbouring villi at $Wo = 0.5$. Instantaneous flow at $t/T = 0.3$ for (A) in-phase ($\Delta\phi = 0$) and (B) out-of-phase ($\Delta\phi = \pi$) villi motion. The colour fields map the magnitude of instantaneous velocity ($|\mathbf{u}|/U_0$). Time-averaged steady streaming flow-fields ($\mathbf{u}^{ss}$) for (C) in-phase ($\Delta\phi = 0$) and (D) out-of-phase ($\Delta\phi = \pi$) villi motion. The colour fields map the magnitude of $|\mathbf{u}^{ss}|/U_0$. The width of the streamlines is proportional to the local velocity magnitude. See also figures S2-S4 in supplementary material.

## 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 for $Wo = 0.5$ in figures 2(C) and (D) for $\Delta\phi = 0$ and $\Delta\phi = \pi$, respectively. The SSF is characterized by counter-rotating vortices just above the villus tips, with amplitudes several orders of magnitude lower than the instantaneous flow velocity.

For the synchronised oscillation case $\Delta\phi = 0$, an outer region of weak unidirectional steady flow develops in the lumen (see figure S4). The strength of this flow increases with increasing $Wo$ in agreement with the Rayleigh streaming theory (Puthumana Melepattu & de Loubens 2022; Tatsuno 1973). No such outer streaming region develops for the case $\Delta\phi = \pi$. The key point here is that for non-propagating contractions, $\Delta\phi = 0$ and $\pi$, irreversible SSF manifests as an inertial phenomenon, and $\mathbf{u}^{ss}(\mathbf{x}) \approx 0$ in the Stokes flow regime when $Wo \ll 1$.

# 5. Pendular-wave contractions and opposing steady streaming

In this section, we investigate the flow dynamics generated by the propagating contractions and relaxations of the intervillous spaces along $-\hat{z}$, as illustrated in Figure 1(D). These correspond to the cases simulated with a phase lags in the range $0 < \Delta\phi < \pi$, for different $Wo$.

## 5.1. Instantaneous flow

Figures 3 and 4 show the flow-fields for $Wo = 0.16$ and $Wo = 2.82$, respectively, for two values of $\Delta\phi$ in panels (A) and (B), at a given time instance. From the figures, we notice the emergence of a boundary layer above the villi tips (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 boundary layer* (MBL) and its associated height is noted as $\ell$. Above the MBL, until the central axis of the lumen at $r = R$, we observe near-uniform unidirectional flow moving along $+\hat{z}$. We term this region as the *advected layer*.

Curiously, the vortical structures in the mixing layer are seen to propagate along $-\hat{z}$,
<figure><img src="image_4.png" /><figcaption>Figure 3: Snapshots of the instantaneous flow field for $\bar{a} = 0.2$ and $Wo = 0.16$, at $t/T = 0.45$, 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 advected layer. Note the absence of the advected layer in (B). See supplementary movies 1 and 2.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4: Snapshots of the instantaneous flow-field for $\bar{a} = 0.2$ and $Wo = 2.82$, at $t/T = 0.45$, 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 advected layer. See supplementary movies 3 and 4.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: 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>

carried by the imposed propagating pendular-wave along the wall, while the axial flow in the advected layer goes opposite, along $+\hat{z}$ (see supplementary movies 1–4). For the case $\Delta\phi = \pi/4$ and $Wo = 0.16$ (figure 3(B)), the mixing layer extends all the way up to $r = R$, and the advected 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 in $Wo$, as can be observed when comparing figures 4 with 3. The decrease of $\ell$ with $Wo$ is shown in supplementary movie 5. For example, in the aforementioned case (at $\Delta\phi = \pi/4$), 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$.

The effect of propagating wave along the villi-wall is felt throughout the lumen zone, irrespective of the MBL height. Figure 5 shows the measured radial velocity $u_r$ at $r = 2H$ in the lumen, in panel (A) $Wo = 0.5$ and (B) $Wo = 2.82$, at increasing 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_z$ (data not shown). This behaviour of the instantaneous velocity is seen both in the mixing and the advected layer, and the wave amplitude is seen to monotonically decreases with increasing radial distance $r$.

Whereas the propagating pendular-wave induces a travelling wave in the lumen (figure 5), the emergence of an advected layer that axially pumps fluid in the opposite direction is an intriguing phenomenon (figures 3 and 4). This behavior is counter-intuitive in the context of momentum transport via travelling waves, and contrasts with canonical peristaltic flow theory (Jaffrin & Shapiro 1971), in which the wave transports the fluid in its direction of propagation.

## 5.2. Steady streaming flow

The axial fluid pumping effect is not limited to the advected layer, but occurs throughout the lumen including in the mixing layer. This is demonstrated by the figure 6 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 $\Delta\phi$ (column-wise). The flow is shown over two neighbour villi (at their mean axial positions) in each simulated system, where the color-field plots the axial component ($u_z^{ss}$) of the SSF velocity. Note that the SSF has a periodicity of $P$ in the axial direction. 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 $+\hat{z}$), descending down through the intervillous spaces and curving above the villi tips. 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.

It is worth noting that the MBL height $\ell$ does not appear to leave any signature on the SSF flow field, as evidenced by the dashed lines in figure 6. We therefore cannot directly correlate the radial variations of the SSF patterns with the emergence of the MBL discussed in section 5.1.

# 6. Mixing boundary layer evolution

In the previous section, simulations revealed a separation of the flow field into a mixing boundary layer above the villi and an advective layer in the center of the lumen. The emergence of an viscous boundary layer, normal to planar boundaries, is a classical phenomenon in pulsating flows at finite inertia (Schlichting 1960; Loudon & Tordesillas 1998). This is typically referred to as the Stokes boundary layer, and scales as $\delta_{St} \sim \sqrt{\mu/(\omega\rho)}$. In this section, we show that patterning the channel with villi leads to departure from the classical theory.
<figure><img src="image_7.png" /><figcaption>Figure 6: 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>

## 6.1. *Mixing boundary layer height*

The MBL height $\ell$ separates the vortical flow region from a unidirectional flow region (advected layer) above it. In our simulations, the axial velocity in the advected layer is always positive, i.e. $u_z > 0$. We use this fact to characterize the height of MBL, by first decomposing
<figure><img src="image_8.png" /><figcaption>Figure 7: Plot of the mixing layer height $\ell$ against the Womersley number $Wo = W/\delta_{St}$, for various $\Delta\phi$, for (A) $\tilde{a} = 0.2$ and (B) $\tilde{a} = 0.1$. Note the plateauing of the curves when $Wo \ll 1$, and their rapid decrease with increasing $Wo$, when $Wo \ge 1$, indicating a flow regime transition. The measured mixing layer heights $\ell$ 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 $\ell$ 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: $\ell^c = 6.66e^{-3\tilde{a}/2}W/W_o$, 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 $\tilde{a}$.</figcaption></figure>

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^{-} \qquad (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, \qquad (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 7 for various $\Delta\phi$, in panels (A) for $\tilde{a} = 0.2$ and (B) for $\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 extends up to the channel centerline ($r = R$), without the emergence of an advected layer.

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$, which depends solely on $\Delta\phi$ and $\tilde{a}$. This is a departure from the classical oscillating boundary layer theory (Schlichting 1960), which predicts that at low $Wo$ (i.e. in the Stokes flow regime), the viscous dissipation layer should encompass the entire channel. 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$. Comparing figures 7(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 markers), the Stokes to inertial transition of $\ell$ occurs at a smaller $Wo$, as compared to that at larger values e.g. $\Delta\phi = \pi/2$ (triangle markers). 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 advected 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 counterclockwise 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 6, panel: $\Delta\phi = \pi/5, Wo = 5.0$) and the vortical instantaneous 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 7(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,
<figure><img src="image_9.png" /><figcaption>Figure 8: 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 $t/T = 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 advected layer. The dashed (magenta) separation line between the mixing and advected layers is reproduced from figure 3(A).</figcaption></figure>

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^{-3\tilde{a}/2}\delta_{St})$. Physically, the $y$ variable represents the MBL height normalized by the viscous penetration layer above the villi, after correcting for spatial attenuation associated with the oscillatory amplitude. The $x$ variable is $\propto 1/Wo$ and quantifies the relative reduction of oscillatory inertial effects, rescaled by the leading-order dependence of the MBL height on the phase lag. For large values of the $x$ variable, the flow is dominated by viscous effects, leading to the power law relation $y \propto x^{-1}$ (blue line), and thus $\ell \propto \Delta\phi^{-2/3}e^{-3\tilde{a}/2}$. Here, the mixing layer height $\ell$ is unaffected by oscillatory fluid inertia and is therefore independent of $\delta_{St}$, identifying the parameter set that falls within the viscous dominated regime.

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^{-3\tilde{a}/2}$. The double power law collapse therefore shows a 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 expressed in terms of $\Delta\phi$ as $\ell^c = 5.5e^{-3\tilde{a}/2}H\Delta\phi^{-2/3}$, and similarly, in terms of $Wo$, as $\ell^c = 6.66e^{-3\tilde{a}/2}W/Wo$. This Stokes-to-inertial critical transition line is overlaid on the unscaled $\ell/H$ data in figures 7(A) and (B) (grey dash-dotted line). The inertial regime is shown by the grey shaded region in all plots in the figure 7.

## 6.2. Competing oscillatory and steady flows

We now look into the physical mechanism that leads to the appearance of distinct mixing and advected layers. Towards this end, we subtract the characteristic flow velocity $U^{AL}$ of the advected layer from the entire flow-field. 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 drdz, \qquad (6.3)$$

making sure that we integrate well above the MBL height, and the radial component is set to $U_r^{AL} = 0$.

Figure 8 shows this shifted velocity field for $\Delta\phi = \pi/2$ and $Wo = 0.16$, with the height
<figure><img src="image_10.png" /><figcaption>Figure 9: Axially averaged $(\cdot)_z$ z-components of (A) the oscillatory axial velocity $\langle u_z^{os} \rangle_z$ and (B) the steady streaming axial velocity $\langle u_z^{ss} \rangle_z$. These measures are computed as functions of radial distance $r/H$, 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$. $\bar{a} = 0.2$. 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>

of the MBL reproduced (dashed magenta line) from figure 3(A). The shifted velocity field consists entirely of counter-rotating, semi-vortical flow structures that now extend up to the central symmetry boundary. This structure is similar to the one seen when the mixing layer extends till the center line (e.g. figure 3(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 axial steady flow and oscillatory components.

The distinct mixing and advected layers seen in the simulations manifest due to the competition between the strength of the oscillatory flow, which generates the counter-rotating vortical structures, and the axial steady flow shown in figure 6. In order to understand the MBL height $\ell$ behaviour, we map the radial dependence of the steady (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)| dt dz, \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 9 compares the strength of these two velocity measurements, for three increasing values of $Wo$ (column-wise panels) and all $\Delta\phi$. In figure 9(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 9(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, the pendular-wave generates a steady 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 so determined by a competition between these two effects. Within the mixing layer the oscillatory axial flow predominates, while in the advected 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 steady axial flow near the channel center, as well as a faster radial decay of the oscillatory axial flow. 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. Origin of irreversible axial flow opposing wave travel

As stated previously, the existence of irreversible flow in low and moderate $Wo$ in a direction opposite to that of the wave travel contrasts with canonical peristaltic flow theory (Jaffrin & Shapiro 1971). In this section we aim to identify the physical mechanism responsible for the steady streaming (irreversible) axial flow opposing the pendular-wave propagation direction, as well as the instantaneous radial pumping. The steady flow component is responsible for the emergent advected layer, as seen in the previous section. To this end, we develop scaling laws that predict the magnitude and direction of the irreversible axial fluxes and radial pumping, and test their validity against numerically integrated measures from the simulation results.

## 7.1. Scaling laws

At first glance, the irreversible axial flow pumping generated by the villi wall studied here is reminiscent of that observed in propagating waves along dense oscillating ciliary arrays (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 to that of net fluid pumping by the ciliary array (Ding et al. 2014). However, individual cilia possess are deformable (an additional degree of freedom) and trace non-reciprocal effective and recovery strokes during their oscillatory cycle. This breaks the time-symmetry of the Stokes equations, and generates irreversible fluid transport. This is clearly not the case for villi, which are rigid structures (Lim et al. 2014) and undergo reversible harmonic oscillations. The non-reciprocal motion of the villi-wall, resulting in steady streaming axial flow even in the Stokes flow regime, is therefore more subtle.

Building on the *Scallop Theorem* (Purcell 1977; Lauga 2011), Najafi (2004) proposed a minimal model for a swimmer in Stokes flow, consisting of three spheres linked in-line, where the inter-sphere gaps undergo a cyclic non-reciprocal sequence of contraction and extension. While the individual spheres follow reciprocal trajectories, the collective behaviour of the system breaks time-symmetry and generates net propulsion in Stokes flow regime (Najafi 2004; Najafi & Golestanian 2005). Referring to figure 1(D), it becomes clear that a similar non-reciprocal sequence is executed by any two adjacent intervillous gaps (and their bounding villi) in our system. The villi wall can thus be viewed as an array of
<figure><img src="image_11.png" /><figcaption>Figure 10: Axial steady-streaming fluxes integrated over (A) the villi-zone $Q_z^{ss,v}$ and (B) over the lumen-zone $Q_z^{ss,l}$, plotted against the phase lag $\Delta\phi$, for increasing $Wo$. Flux quantities are defined as averages taken over the channel length $L_z$. Filled markers (coloured) plot simulations with $\tilde{a} = 0.2$ and empty markers (black) with $\tilde{a} = 0.1$. In both plots, the solid (red) curve is from equation 7.3, and derived from geometric considerations detailed in section A.2. 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>

such non-reciprocally contracting and relaxing gaps, where the additional degree of freedom is introduced via the spatial phase lag $\Delta\phi$. This phase-lagged actuation implies that each downstream (along +z) intervillous gap reaches its minimum or maximum size slightly ahead of its upstream neighbour. This sequence is never reversed in time, ensuring non-reciprocal boundary kinematics.

The contraction-expansion of the intervillous gaps pumps fluid radially, as illustrated by the arrows in figure 1(B). The non-reciprocal motion of adjacent intervillous gaps facilitates the net transfer of some parcels of this fluid from an upstream gap to its downstream neighbour. The villi-wall boundary therefore irreversibly pumps fluid in the axial direction. Based on this mechanistic description, we develop a scaling laws to describe the emergence of the oscillatory radial flow. The rate of expansion and contraction of the gap between two adjacent villi, denoted $\zeta_i$, along the $\hat{z}$ direction is given by,

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

Integrating this equation allows us to scale the peak oscillatory radial flux per intervillous gap, denoted $Q_r^{os,max}$ (see Appendix A.1 for details), which leads to,

$$ Q_r^{os,max} = \frac{2a\omega H}{\pi} \sin\left(\frac{\Delta\phi}{2}\right), \qquad (7.2) $$

Along the axial direction, the fluid trapped in the intervillous gaps undergoes non-reciprocal motion. Taking into account the coupled effect of the gap size and the fluid velocity within the gap, we get the direction and magnitude for the steady streaming axial flux averaged over the domain length $L_z$ (see Appendix A.2). The flux within the villi zone is denoted $Q_z^{ss,v}$ and that within the lumen zone is denoted $Q_z^{ss,l}$, are scaled as,

$$ Q_z^{ss,v} = \frac{a^2\omega H}{2P} \sin(\Delta\phi) \quad \text{and} \quad Q_z^{ss,l} = \frac{a^2\omega(R-H)}{2P} \sin(\Delta\phi). \qquad (7.3) $$
## 7.2. Axial pumping

To confirm our mechanistic scaling arguments, we numerically measure the irreversible axial flux, $Q_z^{ss}$, separately in the villi zone ($0 \le r \le H$) and in the lumen zone ($H < r \le R$) averaged over the channel length $L_z$ as,

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

and

$$Q_z^{ss,l} = \frac{1}{L_z} \int_H^R \int_0^{L_z} u_z^{ss}(z, r) dz dr. \quad (7.5)$$

Figures 10(A) and (B) show the domain averaged axial steady streaming fluxes in the villi and lumen zones, respectively, for various $Wo$. For the irreversible SSF flux in figure 10(A), we see a remarkable data collapse on the scaling equation (7.3) (solid red curve). Here, the data begins to show deviations due to increasing inertial effects only for high $Wo = 5.0$ for $\tilde{a} = 0.1$. The measured data for $Q_z^{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$.

The irreversible luminal flux is plotted in figure 10(B), rescaled appropriately by $\tilde{a}U_0(R-H)$, as a function of $\Delta\phi$. In the viscous-dominated flow regime at low $Wo \ll 1$, the steady-state axial flux data, $Q_z^{ss,l}$, collapses onto the same functional form (solid red curve) given by (7.3). The data-points here for $Wo \le 0.28$ fall slightly below the scaling equation (solid line). Thus, in the Stokes flow regime, we can approximate the total steady streaming axial flux across any height $h$ of the channel as $Q_z^{ss} = Q_z^{ss,v} + Q_z^{ss,l} \approx (Wh \sin(\Delta\phi))/(2PH)$.

However with an increase in oscillatory inertia, the flux data no longer shows a functional collapse. The flux here shows a non-monotonic behaviour with increasing $Wo$ in figure 10(B). Note that $Q_z^{ss,l}$ at $Wo = 5.0$ is lower than that at $Wo = 2.82$, for all $\Delta\phi$. Alt increasing $Wo$ increases the axial velocity of the fluid in the lumen, the increased inertia also limits the radial height to pump this fluid due to dynamic radial confinement of the flow near the villi tips (see figure 6). The flux therefore reaches a maximum at $Wo \approx 2.8$ and falls beyond $Wo > 2.82$.

## 7.3. 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 recirculates, forming counter-rotating semi-vortical flow structures, as seen in figures 3 and 4.

Radial pumping is quantified using radial oscillatory flux. Due to radial flow confinement, mass conservation ensures that a simple integration of the velocity field along the z-axis yields zero net flux. We therefore define the oscillatory radial flux per intervillous gap 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. \quad (7.6)$$

Here, $N = 2\pi/\Delta\phi$ is the number of intervillous gaps along axial periodic length $L_z = NP$. Subtracting $u_r^{ss}(z, r)$ from the instantaneous velocity ensures that we capture the periodically oscillating component of the flux. Unlike in the axial case, the radial flux per intervillous gap, $Q_r^{os}$, is defined to retain its radial dependence and to have equal magnitude in $\pm r$ directions.

Figure 11(A) shows the non-dimensional oscillatory radial flux $Q_r^{os}(r)$ on the x-axis, as
<figure><img src="image_12.png" /><figcaption>Figure 11: (A) The oscillatory radial flux $Q_r^{os}$ per intervillous gap evaluated from (7.6) is plotted (on x-axis) against radial distance $r/H$ (on y-axis) for $\tilde{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$. (B) The peak value of radial oscillatory flux ($Q_r^{os,max}$) per intervillous gap plotted against the phase lag $\Delta\phi$, for increasing $Wo$. Filled markers (coloured) are for $\tilde{a} = 0.2$ and the empty markers (black) for $\tilde{a} = 0.1$. The peak flux data-points collapse solely as a function of $\Delta\phi$, when appropriately rescaled. The simple scaling function is shown with dashed (green) line.</figcaption></figure>

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). The radial flux rises to a peak value just outside the villi tips and then falls to zero at large $r$ for all cases. Unsurprisingly, we observe that the maximum radial flux is obtained for $\Delta\phi = \pi$, corresponding to the case where the intervillous spaces experience the largest strain for a given villus displacement $a$. As $\Delta\phi$ decreases, the peak radial flux value also decreases. Moreover, the radial distance over which the flux decays to zero increases with decreasing $\Delta\phi$; the flux decays most rapidly for $\Delta\phi = \pi$ and most gradually 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_r^{os}$ and its radial dependence is solely controlled by the geometric motion of the villi-wall, set by the phase lag ($\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_r^{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 flux, 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}(r))$. Figure 11(B) plots the rescaled oscillatory radial flux as a function of $\Delta\phi$, for increasing $Wo$. Here, we see a remarkable data collapse onto the scaling form given by equation (7.2) (dashed green line), for all $Wo$ and $\tilde{a}$. Minor deviations are seen for $\Delta\phi = 0$,
where the flux does not remain identically zero for increasing $Wo$. When $\Delta\phi = 0$, the radial oscillatory flux is generated solely through an inertial mechanism, and increases linearly with $Wo$ (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 $Wo$. The peak oscillatory radial flux therefore remains unaffected by fluid inertia ($Wo$) for $\Delta\phi > 0$.

## 8. Effective velocity boundary condition

In this section, we explore an effective velocity boundary condition that captures the primary effects of the pendular-wave along the villi-wall. This boundary condition allows for larger scale simulations of intestinal flow, without explicitly resolving the individual villi. The boundary condition also enables further analytical work in channel flows driven by villi-patterned walls.

Towards this end, we carry out a Fourier decomposition of the instantaneous velocity signal measured at $u(t, r \approx H)$ in our simulations, just above the villi tips, identifying the significant modes appearing at relevant length scales. We therefore propose the following component-wise non-dimensional effective velocity boundary model, $\tilde{u}^e$, for the system, $\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 (8.1) $$

$$ \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 (8.2) $$

where $c_{0-2}$ and $d_{0-2}$, are model coefficients that will depend on $\Delta\phi$ and $\tilde{a}$, as well on the geometric ratios of the domain.

The effective boundary velocity model consists of three terms for both the axial and radial components. 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.

The coefficients in (8.1) and (8.2) are shown in figure S6 as functions of $\Delta\phi$ for $Wo = 0.16$. We obtain the coefficients $c_0$ and $d_0$ by averaging the instantaneous velocity signal measured at $u(t, r \approx H)$ from the simulations over both time and space, for different values of $\Delta\phi$ and $\tilde{a}$. We see that $c_0$ is non-zero, which is necessary for generating irreversible unidirectional axial flow, manifesting the advected layer. The measured $c_0$ coefficient is $\approx 90\%$ of $W/(2P)$, the analytically obtained scaling for velocity within the villi zone in (7.3). For the radial component, $d_0 \approx 0$ since the fluid is radially confined by walls and no net radial flow can occur. The model coefficients can thus be reduced by setting $c_0 = \frac{\tilde{a}W}{2P} \sin(\Delta\phi)$ and $d_0 = 0$.
<figure><img src="image_13.png" /><figcaption>Figure 12: 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 8.1 (solid, red), for four increasing time-fractions. A comparison of the radial velocity $u_r/U_0$ measured at the villus tips with the effective boundary condition is provided in figure S5. The coefficients $c_{0-2}$ for the effective velocity curves plotted here are obtained from the heuristic approximations for the data shown in figure S6.</figcaption></figure>

The rest of the fit coefficients ($c_{1,2}$ and $d_{1,2}$) are estimated by least-squares fitting of the model equations to the velocity signal, for all time values, and subsequently taking their averages. These show a complex dependence on $\Delta\phi$ and $\tilde{a}$ (see figures S6 and S7). Since this dependence is non-trivial, we introduce heuristic approximations for these coefficients to achieve reasonable fits for smaller $\Delta\phi \to 0$. Note that the variation of $c_2$ and $d_2$ coefficients with $\tilde{a}$, introducing villi scale perturbations to the flow, is consistent with the scaling of the MBL height $\ell$ in section 6.1.

We compare the modelled effective axial boundary velocity ($\tilde{u}_z^e$) with the simulated instantaneous axial velocity ($\tilde{u}_z$) at the villi tips ($r = 1.02H$) for $Wo = 0.16$ for two $\Delta\phi$ values, in figure 12. Similar comparison for the radial component given in (8.2) is provided in the supplementary figure S5. Here we have used the functional approximations of the coefficients shown in figure S6 (solid line) to obtain the effective velocity curves in figure 12 and figure S5. The comparisons show excellent agreement with the simulations. We also note that the model fit improves significantly for the smaller value of $\Delta\phi = \pi/4$. We find that the proposed effective velocity boundary condition works well in the Stokes flow regime, but its accuracy decreases for $Wo \gtrsim 1.0$.

It is worth noting that the effective velocity model in (8.1) and (8.2) 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 3 and 4. Furthermore, one can transform (8.1) and (8.2) 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.
## 9. Conclusions

We investigate the time-periodic flow induced by pendular-wave activity in a channel patterned with villi-like micro-structures, using two-dimensional lattice Boltzmann simulations.

We reveal the existence of oscillatory boundary layer that divides the flow into mixing layer nearer the villi tips and an advected layer nearer the channel center. The mixing layer is dominated by asymmetrical semi-vortical flow structures originating at the nodes and anti-nodes of the pendular-wave, and the advected layer shows uniform unidirectional axial flow. The height of this layer scales as $\ell \propto \Delta\phi^{-2/3}$ in the viscous dominated flow regime, and scales as $\ell \propto \Delta\phi^{-1/3}$ in the inertia dominated regime, where $\Delta\phi$ is the phase lag of the villi-wall boundary. We delineate the two flow regimes by achieving a double power law scaling for the measured $\ell$.

We show that the emergence of the advected layer above the mixing layer results from a competition between the decay of the oscillatory axial flow and the unidirectional steady axial flow along the channel. We map the mixing layer decay with increasing oscillatory inertia for the villi-patterned wall as $\delta_{St}^{-1/2}$, a departure from the classical scaling of $\delta_{St}^{-1}$ (Schlichting 1960) above oscillating planar walls. With increased inertia, the axial fluid gets dynamically confined to nearer the villi-tips while the luminal flow velocities increase. This results in a non-monotonic axial pumping effect, with a maximum seen at $Wo \approx 2.8$.

We reveal that the advective layer transports fluid in a direction opposite to that of wave propagation; a counter-intuitive phenomenon that persists even at low Womersley numbers. We demonstrate that this reversal stems from a subtle mechanism rooted in the *scallop theorem* for Stokes flow: the phase lag between adjacent villi imposes non-reciprocal motion on the fluid entrapped in the intervillous spaces, thus breaking the time-reversibility of Stokes equations. By analogy with the simple 2D Stokes swimmer of Najafi (2004), we formulate mechanistic scaling laws that capture the direction and magnitude of the irreversible flux through the channel.

Finally, we develop an effective boundary condition at the villus tips that incorporates both the larger villi-array length scale as well as individual villi scale. The effective condition successfully models both radial and axial villi-tip velocities for $Wo \le 1$, enabling coarse-grained modelling of intestinal flows at the organ scale without the need to explicitly resolve individual villi.

These results shed new light on the flow consequences of both propagating longitudinal contractions and the presence of villi in the small intestine (Melville *et al.* 1975; Lammers 2005; Lentle *et al.* 2012; de Loubens *et al.* 2013; Fullard *et al.* 2014). On the one hand, villi have long been assumed to play a passive role in absorption, primarily by increasing the available surface area (Strocchi & Levitt 1993). On the other hand, seminal theoretical modelling studies have demonstrated that efficient mixing in the small intestine requires conditions akin to a longitudinal array of perfectly mixed stirred-tank reactors, with slow axial transport of reactants between them (Penry & Jumars 1986, 1987). While it is well established that pendular contractions promote local dispersion via shear dispersion mechanisms (de Loubens *et al.* 2013; Fullard *et al.* 2014), our results highlight that villi themselves could play an active role in enhancing transport and mixing. However, our study also reveals that the physiological basis of intestinal motility at the scale of villi remains largely undocumented. In particular, quantitative data describing the amplitude, frequency, and spatio-temporal organization of villus-scale motions are still lacking. Identifying these missing physiological inputs is essential to advance our mechanistic understanding of transport and mixing phenomena in the intestinal lumen. Our simulations also provides a foundation for future investigations into physiologically relevant effects, such as the influence
of non-Newtonian and heterogeneous digesta, as well as radially asymmetric contractile patterns, on transport mechanics in the gastrointestinal tract.

Future work should extend these 2D simulations to 3D geometries to quantify how villus-scale architecture influences intestinal transport. The transport mechanisms identified here are expected to hold in 3D, as the scaling law is based on inter-villus mass conservation, independent of villus shape or dimensionality. When moving to a cylindrical geometry, we expect from mass conservation arguments that circumferential confinement will enhance the advective layer near the tube center. Among geometric parameters, the lateral spacing between villi—both longitudinal and circumferential—appears particularly critical, controlling the thickness of the steady streaming boundary layer (Puthumana Melepattu & de Loubens 2022), and should therefore guide future 3D studies. These extensions will enable a more complete quantitative understanding of mixing and axial transport in realistic intestinal geometries.

Beyond these physiological considerations, these findings suggest that bio-inspired “pendular-wave” villus-wall motion could be harnessed for microfluidic applications. Much like artificial cilia arrays (Shields et al. 2010; den Toonder et al. 2008), such motion could enable efficient flow control in microfluidic systems. In this context, varying the villi geometry or spacing ($R/H$, $P/W$) could provide an additional means to tune the fluid fluxes, as seen in related studies (Puthumana Melepattu & de Loubens 2022; Hall & Clarke 2020). This opens new avenues for biomimetic design, where pulsatile boundary and pressure forcing offer an additional degree of freedom for flow control.

**Supplementary data.** Supplementary movies are available. The simulation code is open source and available at Vernekar et al. (2025).

**Acknowledgements.** RV thanks Irina Ginzburg for her advice 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 projects TransportGut, ANR-21-CE45-0015 and TABAG, ANR-20-CE30-0001-01.

**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;

**Author contributions.** RV and FA performed numerical computations. RV and MG developed theoretical scalings. DIYM and ST provided experimental inputs. RV and CdL carried out analysis of the results. RV authored the computational code. CL, ST and CdL supervised the research and managed funding. All authors discussed the results and contributed to the final manuscript.

## 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 (A3) $$

### A.1. 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^{\text{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 (A4) $$

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^{\text{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 (A5) $$

which gives the non-dimensional oscillatory radial flux as,

$$ \tilde{Q}_r^{\text{os,max}} = \frac{2}{\pi} \sin\left(\frac{\Delta\phi}{2}\right), \qquad (A6) $$

### A.2. 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) \quad (A7) $$

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 (A8) $$

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

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

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} \bar{M}_z = \frac{a^2 \omega H}{2P} \sin(\Delta\phi). \quad (A10) $$

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 (A11) $$
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.

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.

CREFF, JUSTINE, MALAQUIN, LAURENT & BESSON, ARNAUD 2021 In vitro models of intestinal epithelium: Toward bioengineered systems. *Journal of Tissue Engineering*.

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.

FULLARD, LUKE, LAMMERS, WILLEM, WAKE, GRAEME C. & FERRUA, MARIA J. 2014 Propagating longitudinal contractions in the ileum of the rabbit – Efficiency of advective mixing. *Food & Function* **5** (11), 2731–2742.

GARIC, MARTIN, VERNEKAR, ROHAN, MARTÍN, DÁCIL I YÁNEZ, TANGUY, STÉPHANE, DE LOUBENS, CLÉMENT & LOVERDO, CLAUDE 2025 Intestinal villi and crypts density maximizing nutrient absorption. *arXiv preprint arXiv:2507.03472*.

GINZBURG, IRINA 2025 The Lattice Boltzmann Method with Deformable Boundary for Colonic Flow Due to Segmental Circular Contractions. *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.

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.

HARDACRE, ALLAN K, LENTLE, ROGER G, YAP, SIA-YEN & MONRO, JOHN A 2018 Predicting the viscosity of digesta from the physical characteristics of particle suspensions using existing rheological models. *Journal of The Royal Society Interface* **15** (142), 20180092.

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, 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, 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.

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.

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.

PENRY, DEBORAH L & JUMARS, PETER A 1986 Chemical reactor analysis and optimal digestion. *Bioscience* **36** (5), 310–315.

PENRY, DEBORAH L & JUMARS, PETER A 1987 Modeling animal guts as chemical reactors. *The American Naturalist* **129** (1), 69–96.

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.

VERNEKAR, ROHAN, GARIC, MARTIN, DE LOUBENS, CLÉMENT, LOVERDO, CLAUDE & TANGUY, STÉPHANE 2025 3dintestinalflow: Computational code for simulating intestinal flow with villi-induced motility. HAL archive, hal-05155129, version 1.

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.

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.

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