Shear-induced dispersion in peristaltic flow

Brownian tracers or solutes in a quiescent fluid disperse as a result of molecular diffusion only. Unidirectional flows such as pipe flows, however, stretch and enhance tracer concentration gradients, leading to an increased effective diffusivity at long times. In a landmark paper, Taylor¹ built upon this basic picture to analyze the dispersion of a solute in a tube. He arrived at a cross-sectionally averaged advection–diffusion equation for the tracer concentration with a shear-dependent diffusivity known as the dispersivity. This effective dispersion coefficient was found to scale as Pe², where Pe = Ua/D is the Péclet number of the flow expressed in terms of the mean velocity U, channel radius a, and molecular diffusivity D. Taylor’s calculations were based on strong intuition that was further formalized by Aris² in the so-called method of moments. In the context of transport in a tube, the Taylor–Aris prediction is an asymptotic result and is only valid once the tracers have had sufficient time to sample all the transverse positions in the channel cross section. Since these seminal models, many other methods have been proposed to calculate the effective dispersivity, which involve asymptotic expansions,³ Frankel and Brenner’s generalized Taylor dispersion (GTD) theory,⁴ center-manifold reduction,^5,6 and, more recently, a formulation that combines the method of moments with Dirac’s bra–ket formalism.⁷ All of these methods rely on a well-known asymptotic technique consisting of eliminating fast modes in a problem to arrive at a simpler equation for the long-time behavior of a slow mode.⁶ 

The concept of a simplified cross-sectionally averaged transport equation is appealing since in many applications one is primarily interested in determining asymptotic transport properties such as the mean velocity or mean square displacement. As a result, the classical analysis has been extended and applied to various problems involving dispersion in turbulent flows,⁸ solute transport in pipes with reactive walls,⁹ cross-flow problems to model sedimenting particles or filtration,^10–12 geophysical flows,¹³ and microfluidics.^14–16 

A particularly relevant extension of Taylor dispersion in the present context applies to transport through porous media. The analysis of solute transport in porous media flows was formalized by Brenner,¹⁷ who proposed a general theory for dispersion in a spatially periodic matrix based on GTD, which we briefly review in Sec. II B. This GTD theory can be used to calculate the effective dispersivity by the solution of a conservation equation and quadratures over the representative unit cell defining the periodic lattice. This technique has since been used extensively to study shear dispersion in porous materials,^18,19 electrophoretic²⁰ and pressure-driven²¹ flows in periodic and serpentine channels,²² and periodic networks²³ in the context of microfluidic applications. These various models come under the purview of “macrotransport theory,” which aims to derive asymptotic equations for measurable long-term quantities from the governing equations of microscopic field variables.

In this work, we make use of Brenner’s GTD theory for porous media¹⁷ to study shear-enhanced dispersion under peristalsis. The flow is driven by a prescribed periodic wave train on the flexible walls of a two-dimensional channel, resulting in net unidirectional pumping. Peristaltic pumping is quite ubiquitous in biological processes, and examples include transport in the ureter, in the digestive tract, and in certain types of blood vessels. The mechanism has also been exploited to pump fluid in various microfluidic devices.^24–27 Understanding the transport of passive Brownian solutes in these flows is, therefore, a problem of fundamental interest. There has been a number of numerical investigations of transport in perstaltic flows involving particle simulations and a dynamical systems approach.^26,28 In contrast, theoretical investigations have been scarce²⁹ and have focused on the limit of long wavelength deformations.

In the absence of flow, the transport of Brownian tracers in periodic geometries or periodic potential landscapes has been analyzed in the context of Fick–Jacobs (FJ) theory.³⁰ This approach models effective diffusive dynamics of particle positions in entropic potentials.³¹ Macrotransport GTD theory has been proposed as an alternative to FJ theory to understand force-driven transport through entropic barriers.^32,33 The theory and numerics developed in the present paper capture this classical limit of pure diffusion as a special case but also extend it to account for peristaltic flow.

The paper is organized as follows: we start by laying out the problem definition and governing equations and then review the basics of GTD theory in spatially periodic porous media in Sec. II. We solve the flow and dispersion problems in the long-wave limit using the lubrication approximation in Sec. III, where we obtain analytical expressions for the dispersivity in various asymptotic limits. We then generalize the results to arbitrary channel aspect ratios in Sec. IV, focusing on the Stokes limit where we use the boundary integral method for flow calculations. Solving for the dispersivity in that case involves the solution of a conservation law using a finite-volume method, details of which are also presented in Sec. IV. We summarize the results and discuss possible extensions in Sec. V.

We analyze the solute transport in an infinite two-dimensional channel whose walls deform periodically in the transverse direction according to a sinusoidal traveling wave [Fig. 1(a)]. We denote variables in the fixed laboratory frame with tildes $̃$⁠, with the position vector defined as $R̃=(X̃,Ỹ)$ with respect to a fixed origin located on the channel centerline. The half-width of the channel from the centerline is given by

where 2h₀ is the mean channel width, c is the wave speed, λ is the wavelength, and γ is a geometric parameter controlling the amplitude of deformation. The location of the two channel walls is then given by $Ỹ=±h̃(X̃,t̃)$⁠. We also introduce A = 2h₀λ as the area of the channel enclosed by one wavelength. In this fixed frame, material points on the walls move periodically in the $Ỹ$ direction with velocity $±∂h̃(X̃,t̃)/∂t̃$⁠.

It is convenient to perform a Galilean transformation to a moving frame that translates with the wave speed c in the $X̃$ direction.³⁴ To this end, we define the new coordinates

where $(x^,y^)$ are unit vectors in the streamwise and transverse directions. In the moving frame, the geometry of the channel is fixed with respect to time,

and a velocity of − c exists at the walls in the X direction [Fig. 1(b)]. We also note that the partial derivative with respect to time is transformed according to

Denoting the fluid velocity in the fixed reference frame as $ũ=(ũ,ṽ)$⁠, the corresponding velocity u = (u, v) in the moving frame is simply given by

and is independent of time. It satisfies the steady incompressible Navier–Stokes equations,

Periodic boundary conditions apply in the X direction. At the walls, the no-slip condition in the fixed frame requires u(X, Y) to have a slip velocity of –c,

A boundary condition on v is obtained from the no-penetration condition $n^⋅u=0$⁠, where $n^$ is the unit normal. At the top wall,

Combined with Eq. (7), this provides the condition

Equations (6), along with boundary conditions (7) and (9), entirely specify the flow problem. The formulation presented here is classic, going back to the seminal work of Shapiro et al.³⁵ Two different solution approaches are presented in the following, based on a long wavelength approximation in Sec. III and on the boundary integral method at zero Reynolds number in Sec. IV. Having determined the flow field, a quantity of interest is the net flow rate or pumping rate. In the fixed laboratory frame,

where

denotes the flow rate in the moving frame and is a constant independent of both position and time.³⁶ 

We analyze the transport of a single Brownian tracer with instantaneous position $R̃(t̃)$ in the fixed laboratory frame. The tracer is advected by the fluid flow and also diffuses with molecular diffusivity D. The statistics of $R̃(t̃)$ are governed by the probability density function $ψ̃(R̃,t̃)$⁠, which satisfies the Fokker–Planck equation,

where the flux $J̃$ is given by

The right-hand side of Eq. (12) captures the initial condition, with the tracer located at $R̃=R̃′$ at $t̃=0$⁠. The Fokker–Planck equation is subject to the no-flux condition on the walls of the channel, $n^⋅J̃=0$⁠, and to the normalization condition,

where the integral is over the entire area Ω_∞ of the channel.

Our goal in this paper is to predict the long-time behavior of the ensemble-averaged tracer position and variance, respectively, defined as

and

At long times (⁠$t̃≫h02/D$ and $t̃≫λ2/D$⁠), both of these moments are expected to grow linearly with $t̃$⁠,

The vector $Ũ*=Ũ*x^$ characterizes the mean transport velocity, whereas the second-order tensor $D̃*=D̃*x^x^$ describes the mean dispersivity. By symmetry, they each only involve one scalar quantity: the mean velocity $Ũ*$ and dispersion coefficient $D̃*$⁠, respectively.

As for the flow problem, we move to the reference frame translating with the peristaltic wave by applying the transformation of Eqs. (2) and (5). By virtue of Eq. (4), we find that the Fokker–Planck equation (12) in the moving frame remains unchanged,

where ψ(R, t) is now the probability density function of finding the particle at position $R=R̃−ct̃x^$ in the new frame, and the form of the flux also remains identical,

Equations (15) and (16) for the ensemble-averaged position and variance transform as

and, therefore, the mean velocities and dispersion coefficients in the two frames are related by

We now proceed to outline the solution for U^* and D^* in the moving frame based on GTD theory.³⁷ More precisely, we apply Brenner’s macrotransport scheme for spatially periodic media,¹⁷ which we specialize to the singly periodic geometry of the channel along the X direction. The model also closely follows the calculations of Yariv and Dorfman²⁰ for electrophoretic transport in channels with periodically varying cross sections. We simply quote the relevant results here without proof, and the reader is referred to Brenner’s seminal paper¹⁷ for more complete derivations.

We recall that in the moving reference frame, the channel geometry is fixed and given by Eq. (3) and is periodic along the X direction with period λ. We can, thus, think of the channel as being generated by the discrete translation of a unit cell of length λ, and it is convenient to index each cell by an integer $n∈Z$⁠, with the reference n = 0 corresponding to the unit cell where the origin is located (Fig. 2). The global position R = (X, Y) of a point tracer inside the channel can then be decomposed as $R=nλx^+r$⁠, where n specifies the cell in which the tracer is located and r = (x, y) is a local coordinate inside that cell. In other words, the instantaneous position R of the tracer in the channel is equivalently specified by the pair (n, r), which forms the basis of Brenner’s theory. The probability density function ψ(R, t) is then expressed as ψ(n, r, t). Note also that, by periodicity, the fluid flow in the channel is only a function of the local coordinate, u(R) = u(r).

The first step toward solving for U^* and D^* consists in determining the asymptotic local probability density $ψ0∞(r)$⁠, which describes the long-time probability density function for finding the tracer at local position r, irrespective of its cell index n. We obtain $ψ0∞$ by summing the probability density function ψ(n, r, t) over all values of n,

The governing equation for $ψ0∞(r)$ can be inferred from Eq. (18) to be

with flux

This is subject to the no-flux condition at the walls, periodicity requirements at x = 0, λ, and the normalization condition

where Ω_c denotes the domain of a unit cell. Because the flow is incompressible, the solution to this problem is trivial,

where A = 2h₀λ is the area of the unit cell. The mean transport velocity is then given by

The contribution from the diffusive flux in Eq. (25) can be eliminated using the divergence theorem and boundary conditions. After projection along the x direction, we obtain

In other words, the asymptotic mean velocity is simply given by the cross-sectionally averaged fluid velocity.

Solving for the dispersion coefficient requires determining the so-called B-field, which is the local scalar defined as²⁰ 

It characterizes deviations in streamwise tracer transport from the mean velocity and can be interpreted as a dispersion potential. It satisfies the advection–diffusion equation,

subject to the no-flux condition at the walls, $n^⋅∇rB=0$⁠, and to a jump condition at the periodic boundaries of the unit cell,

Once the B-field is known, the dispersion coefficient is simply calculated as

The solution of the B-field is only determined up to a constant, whose value does not affect the result for D^*.

In this section, we first present an asymptotic solution for the flow and dispersion problems based on the classical lubrication approximation for small aspect ratio channels with h₀/λ = ε ≪ 1. Our calculation for the flow field is classic and closely follows the work of Latham,³⁴ whereas the solution of the dispersion problem shares similarities with the work of Yariv and Dorfman.²⁰ 

Following lubrication theory, we non-dimensionalize the coordinates and flow variables using the following scales:

The dimensionless Navier–Stokes equations read

where the Reynolds number is defined as Re = ρch₀/μ based on the channel half-width. At the top wall located at y = h(x) = 1 + γ sin(2πx), the boundary conditions are

By symmetry, we only need to solve for the flow in the upper half of the channel and, therefore, apply symmetry boundary conditions at the centerline,

The dimensionless flow rate, scaled by h₀c, is given by

From here on, we exclusively use dimensionless variables unless specified otherwise.

We solve the Navier–Stokes equations by seeking a regular perturbation expansion in powers of ε ≪ 1,

We shall later see that, for the purpose of estimating dispersion, only the leading-order axial velocity field is needed. Here, we provide the leading order axial and the transverse [$O(ε)$] velocity fields. At leading order in ε, Eq. (37) shows that the pressure is at most a function of x: p₀ = p₀(x). Since the leading order transverse velocity field is of order $O(ε)$⁠, the x momentum equation (36) simplifies to

Using the no-slip and symmetry boundary conditions, this integrates to the well-known parabolic profile typical of lubrication problems,³⁸ 

The above equation must integrate to the leading-order flow rate Q₀ according to Eq. (40). On simplifying, this leads to an expression for the local pressure gradient,

Given the pressure gradient, we can define the net pressure rise over one unit cell of the channel,³⁴ 

where we introduce the notation

If there is no imposed pressure gradient, then the pressure rise must be identically zero by periodicity. This is a case of pure peristalsis where the flow is entirely driven by the imposed motion of the walls. Dispersion in pressure-driven flow is already well understood since the work of Taylor¹ and Aris,² so we focus on pure peristaltic flow and set PR = 0 in this work. Under this condition, the flow rate is entirely determined from the geometry of the walls,

from which we find the mean transport velocity as

Combining Eqs. (44) and (45) yields an expression for the axial velocity in terms of the flow rate and geometry,

Substituting this expression into the continuity equation (35) and integrating yield the leading-order transverse velocity of order $O(ε)$⁠,

This completes the solution to the leading order hydrodynamic problem.

We apply the same lubrication scalings as in Eq. (34) to the governing equation (31) for the B-field, also scaling the variable B with the wavelength λ. After simplifications, the governing equation becomes

where Pe = cλ/D is the Péclet number and measures the relative strength of advection over diffusion. The no-flux boundary condition at the top wall is expressed as

whereas the symmetry boundary condition at the channel centerline implies

Finally, the jump condition (32) becomes

The axial dispersion in the long wavelength limit is then simply given by the following quadrature over the unit cell:

where D^* has been scaled by the molecular diffusivity D.

We seek a solution for B as a regular perturbation expansion in the small parameter ε²,

It should be noted that the validity of this expansion is more subtle than for the flow problem. Indeed, all the terms appear to be of order $O(ε2)$ on the right-hand side of Eq. (52) due to the choice of non-dimensionalization. However, for large Péclet numbers, the advective terms will no longer be $O(1)$ so that the appropriate condition for the validity of the expansion becomes ε²Pe ≪ 1. Substituting the expansion (57) into Eq. (52), we proceed to solve for B at the first two orders.

Solution at $O(0)$⁠: At the zeroth order, we simply have

with boundary conditions

This implies that B₀ = B₀(x), so to leading order the B-field is independent of the transverse coordinate. In order to solve for this leading order field, we must consider the next order term in the expansion.

Solution at $O(ε2)$⁠: Upon collecting terms of $O(ε2)$ in Eq. (52), the governing equation for B₁ is obtained as

where u₀ is the zeroth order solution of the hydrodynamic problem obtained in Eq. (50). The boundary conditions at this order are given by

We integrate Eq. (60) with respect to y over the width of the channel. Using the boundary conditions and recalling the definition (40) for the flow rate, we obtain a differential equation for B₀,

where we have substituted $U0*=Q0/2$⁠. This is subject to the jump condition (55), which can be recast in the form

Equation (63) can be integrated once to obtain

with integrating factor

The integration constant C is obtained by the application of the modified jump condition (64),

Note that the knowledge of $B0′$(x) in Eq. (65) is sufficient to determine the leading order dispersivity using Eq. (56). After simplifications, we, indeed, find

Equation (68), along with (65), provides an expression for the dispersivity, which must be evaluated by numerical quadratures as we do in Sec. III F. This numerical integration can prove to be a challenge at large values of Pe. For this reason and to gain better insight, we complement Eq. (68) with asymptotic results in the two limits of small and large Péclet numbers.

In the limit of Pe ≪ 1, solving for B₀ involves a regular perturbation problem. We seek an expansion of the form

which can be substituted into Eq. (63). The leading order differential equation is simply

This is an exact differential and can be integrated to obtain

where the jump condition ⟨b′⁽⁰⁾⟩ = −1 is used to find the integration constant K = −⟨1/h⟩⁻¹. Substituting this expression into Eq. (68) yields

This shows that the shape modulation in the channel affects the diffusivity to leading order in the absence of any flow. For a straight channel with γ = 0, the leading order dispersivity is simply unity. However, the diffusivity in the limit Pe → 0 drops below one for any value of γ > 0, and it tends to zero when γ → 1 corresponding to the case of fully closed pores. We discuss this behavior in more detail in Sec. III F.

In this limit, we introduce the small parameter δ = 1/Pe and rewrite Eq. (63) as

At first glance, this may seem like a singular perturbation problem involving a boundary layer of thickness δ.³⁹ However, we will find that the asymptotics are simplified by a regular perturbation expansion in the small parameter δ and that the leading order solutions do satisfy all the boundary conditions. Hence, we seek an expansion of the form

At zeroth order, Eq. (73) simply reads

which automatically satisfies the jump condition ⟨b′⁽⁰⁾⟩ = −1 since ⟨h⟩ = 1. At the next order $O(δ)$⁠, Eq. (73) yields

On using the zeroth order solution, we find

which also satisfies the appropriate jump condition due to the periodicity of h. Substituting Eq. (74) into Eq. (68) then provides the dispersion coefficient. After simplifications,

which, for the sinusoidal geometry of Eq. (3), gives

We recall that the large Pe asymptotics in the long wavelength limit are valid when 1 ≪ Pe ≪ ε⁻².

We use Eq. (65) in conjunction with Eq. (68) to calculate the effective dispersivity for a range of Pe numbers and for several values of the geometric parameter γ. For large Pe numbers, the numerical quadratures in Eq. (65) become challenging, and we resorted to use chebfun⁴⁰ and solved Eq. (63) instead. We validate our calculation against Brownian dynamics simulations following the algorithm described in the  Appendix. Figure 3 shows a typical distribution of tracers in one of these simulations, as well as the evolution of the mean-square displacement used to estimate D^*. As we show next, our theoretical calculations match results from simulations very well, as well as asymptotic predictions in the low and high Péclet number limits.

The dependence of the diffusivity $D0*$ on Pe is illustrated in Fig. 4 for different values of γ. As a baseline to understand the results, let us first consider the case γ = 0, corresponding to a straight channel. In that case, $D0*=1$ for all values of Pe, i.e., the dispersion is of purely molecular origin. This is to be expected since the walls are fixed in that case and there is no flow in the channel. As γ is increased, the dispersivity at low Pe falls below unity. In the limit of Pe → 0, the process is also purely molecular and the decrease in effective diffusivity can be understood using Fick–Jacobs theory:³⁰ when γ > 0, the narrow pores act as “entropic barriers” as fewer transverse positions are available there, whereas the wide regions of the channel act as “entropic traps” where the Brownian motion is effectively slower compared to a straight channel. At low Péclet numbers, this has the effect of hindering axial dispersion, with the long-time dispersivity decreasing with decreasing pore opening (increasing γ).³³ In the limit of γ → 1, the long-time dispersivity tends to zero at low Pe, corresponding to the case of completely closed pores with the tracer particles remaining trapped inside the unit cells where they are released.

The scenario becomes quite different as Pe becomes of order $O(1)$ and higher. A net increase in the dispersivity takes place due to the flow, and the trends with respect to pore opening γ are reversed compared to low Péclet numbers. At large Pe, the transport is, indeed, dominated by the mechanism of Taylor dispersion by which the shear generated by the channel deformations couples with cross-streamline molecular diffusion to enhance axial dispersion. This explains both the increase in $D0*$ with Pe and γ since the characteristic flow shear rate scales with Pe and is enhanced in channels with small pores due to incompressibility. We note the peculiar limit of γ → 1 at high Péclet numbers, where $D0*$ does not tend to zero as expected for closed pores and in fact asymptotes to 1 + 3γ²/2. This can be explained by a divergence of the shear rate inside the pore in this limit.

We have solved the dispersion problem in Sec. III in the long wavelength limit of ε ≪ 1, ε²Re ≪ 1, and ε²Pe ≪ 1. However, biological peristaltic pumping often occurs in channels with finite aspect ratios, where the flow topologies can become more complex and involve recirculation bubbles⁴¹ as we show below. The Reynolds number, however, remains low to moderate in many physiological flows,⁴² with Shapiro³⁵ estimating that it is of order unity in the ureter. With these applications in mind, we now use numerical simulations to analyze the case of arbitrary channel aspect ratios in the Stokes regime, with no restrictions on Pe.

For flows with negligible inertia (Re → 0), the velocity field in the moving frame satisfies the Stokes equations,

Following Pozrikidis,^41,43 we solve these equations using the boundary integral representation, which expresses the fluid velocity at a point r on the channel wall as

where f is the distribution of surface tractions on the channel walls. The tensor G is known as the Stokeslet and is the fundamental solution to the Stokes equations for a point force, and T is the associated Green’s function for the stress. Exploiting the periodicity in the x direction, along with the symmetry of the flow about the centerline, we follow Pozrikidis^41,44 and recast Eq. (81) as an integral over one period of a single wall by making use of the appropriate Green’s functions. Upon applying the no-slip and no-penetration boundary conditions, Eq. (81) provides an integral equation for the tractions f, which, once known, are used to determine the velocity at any point inside the channel.

We discretize the top wall of the channel using N + 1 collocation points and N linear elements. The unknown tractions are assumed to be constant over each element, whereas the velocities are assumed to vary linearly. The discrete form of the boundary integral equation (81) at the center r_m of element m is

The matrix A^nm and vector B^mn are given by

where S_n are the linear elements used to approximate the boundary. The integrations over each of the elements are performed using a Gauss–Legendre quadrature rule with 16 quadrature points. The linear system resulting from the discretization of the boundary integral equation results in a dense matrix that is solved using LU decomposition. All the results shown below were obtained using N = 256 collocation points. Once the tractions are known, the boundary integral equation is used to calculate the velocity at any interior point in the domain.

We now turn to the numerical calculation of the dispersion coefficient, which requires solving Eq. (31) for the B field using a finite volume method (FVM). Given the jump condition of Eq. (32), we first find it convenient to define the new variable B^* = B + x, which is now periodic in the x direction. It satisfies the modified governing equation

where the last term −Pe u arises from the change of variable. The aspect ratio ε is no longer constrained to be small. In order to solve this equation, we use a coordinate transformation that maps the physical domain (half unit cell) to a square by scaling the vertical coordinate by the local channel half-width (see Fig. 5),

The partial derivatives are related through the Jacobian matrix,

After transformation, Eq. (85) is integrated over a square grid cell to yield a system of algebraic equations for B^*. The coordinate transformation results in additional non-orthogonal terms that are handled using an implicit formulation.⁴⁵ 

The modified field B^* is subject to periodic boundary conditions in the ζ direction,

and to the symmetry boundary condition at the centerline,

The boundary condition at the top wall is inferred from the no-flux condition, $n^⋅∇rB=0$⁠, which becomes

Having solved for B^*, we obtain the dispersivity by quadrature,

where the dimensionless area is A = 2ε.

The discretized problem for B^*, as presented here, results in a singular matrix. This should not come as a surprise since B^* is only defined up to a constant, and only its gradient is needed to obtain the dispersivity. To specify B^* uniquely, we must impose a solvability criterion, and following the work of Bolster et al.,²¹ we impose a constraint of zero mean for B^* using a Lagrange multiplier. The no-flux boundary condition [Eq. (90)] on the top wall is satisfied iteratively using the method of deferred correction,⁴⁶ which uses solutions from previous iterations to calculate Φ(ζ) and update the source terms resulting from boundary conditions. In most simulations, convergence to 10⁻⁹ is achieved in ∼5 to 40 iterations for a grid of size 200². In the diffusion dominated regime for channels with large aspect ratios, convergence becomes challenging, and in that case, we use a relaxation method with a relaxation parameter of ∼0.5.

We first validate our numerical method by comparison with the long-wave asymptotic results derived in Sec. III, which are valid in the limit of ε²Pe ≪ 1. This limit is difficult to capture using the FVM method, especially for moderate to large Péclet numbers, which require increasingly small aspect ratio channels. The coordinate transformation in the FVM method makes computations challenging for small aspect ratio channels and small pore openings. On the other hand, as the aspect ratio is decreased, the boundary integral method also requires a higher number of quadrature points, which proves to be time consuming. We first compare FVM results against asymptotics in the diffusive limit in the absence of any flow, and excellent agreement is found, as shown in Table I, for an aspect ratio of ε = 0.05. A comparison up to $Pe∼O(5)$ was also performed (not shown), with agreement within 5%. For the latter, the lubrication solution for the flow field was used to calculate the advective terms in the FVM formulation.

We also validate the method for finite aspect ratio channels by comparison with Brownian dynamics simulations following the algorithm outlined in the  Appendix. In these simulations, the velocity field is computed using the boundary integral method and is pre-tabulated on a grid. Linear interpolation is then used to determine the instantaneous position of the tracers during the solution of the Langevin equation (A1). As shown in Table II, both methods show good agreement over a wide range of Péclet numbers.

We now discuss some aspects of dispersion in finite aspect ratio channels with the help of the numerical techniques outlined above. We first consider the diffusive limit in the absence of any flow (Pe = 0), which reduces to the well-studied problem of diffusion in a corrugated channel.^21,33 In this case, our previous discussions based on Fick–Jacobs theory and the concept of entropic barriers still hold. When Pe = 0, the dispersion is entirely controlled by the channel aspect ratio ε and width modulation γ, and trends in terms of these two parameters are presented in Table III. For a fixed width modulation γ, channels with larger aspect ratios ε have a decreased value of their dispersivity. This trend, which is consistent with predictions of Fick–Jacobs theory,³³ may seem counterintuitive as increasing ε for a fixed γ increases the width of the pores, but it also leads to an increase in the size of the entropic traps between pores. In addition, increasing the width modulation for a given ε also causes a sharp decrease in D^* for the same reason as in Fig. 4: as γ approaches 1, the pore openings become increasingly small and strongly hinder longitudinal transport.

Turning to the effects of flow, we first show in Fig. 6 some typical streamlines for different combinations of γ and ε. We find it convenient to consider streamlines in the moving frame, where they are left-right symmetric and lend themselves to an easier interpretation of our results. We recall that even though the flow is different in the fixed frame, the dispersion is the same in both frames. While, in the lubrication limit, the flow fields are always locally parabolic according to Eq. (50), the flow features become more interesting as the aspect ratio ε and width modulation parameter γ are increased. Indeed, recirculation bubbles start to appear, as shown in Fig. 6, which are reminiscent of those occurring in shear flow over cavities.⁴⁷ These recirculation regions are expected to have an effect on tracer dispersion, as noted in previous studies on pressure-driven Stokes flow.⁴⁸ 

The dependence of the dispersivity on the Péclet number for various channel geometries is shown in Fig. 7. The low-Pe behavior is identical to the diffusive limit discussed previously. The effect of the flow becomes noticeable for Pe ≳ 1 − 10, where the dispersivity starts deviating from the low-Pe asymptote and shows an increase with Pe. This change is due to shear-enhanced dispersion, which starts dominating the effect of molecular diffusivity. Finally, in strong flows, the dispersivity eventually increases as Pe², corresponding to the classical Taylor dispersion limit.^1,2 The behavior of the effective dispersivity with Pe is similar to that observed in past studies in the context of dispersion in porous media.^18,49,50

Figure 7 also illustrates the dependence of D^* on ε and γ. The trends at low Pe are consistent with the results of Table III for the diffusive limit. In stronger flows, however, the effects of ε and γ on D^* are reversed: channels with the largest values of ε and γ, i.e., wide channels with large amplitude modulation exhibit the strongest long-time dispersion. In particular, we find that the increase in dispersivity with ε and γ at high Péclet numbers persists after the onset of recirculation bubbles in the flow fields (Fig. 6). This may come as a surprise initially as one may expect the recirculation bubble to act like a hydrodynamic trap. However, the increase in dispersivity can be appreciated from a Lagrangian point of view. Indeed, flows with recirculation bubbles exhibit stronger cross-streamline velocity gradients, and trapping inside recirculation zones increases the residence time of tracers inside the unit cell. Both of these effects, in the presence of molecular diffusion, act to, further, enhance smearing of Lagrangian particle clouds. As a result, the asymptotic dispersivity is increased even though the pre-asymptotic time required to reach the diffusive limit is longer due to the longer residence time inside the unit cell.²¹ 

We have presented a model based on generalized Taylor dispersion theory for the long-time dispersion of a passive tracer transported in peristaltic flow resulting from periodic surface undulations of an infinite two-dimensional channel. The model was first analyzed theoretically in the long wavelength limit, followed by numerical calculations using a boundary integral formulation coupled to a finite volume method in the case of finite aspect ratio channels in the Stokes regime. Our results were also validated against Brownian dynamics simulations and were shown to share many similarities with past studies on dispersion in rigid corrugated channels. In the diffusion dominated regime (low Péclet number limit), the small pores created by channel contractions act as entropic barriers, while the channel bulges act as entropic traps, leading to a reduction in the long-time axial dispersivity in comparison with a straight channel. When a flow is applied, the dispersivity is increased by a mechanism analog to classical Taylor–Aris dispersion, with an $O(Pe2)$ dependence of the dispersivity in strong flows. Interestingly, the trends in geometry are opposite to those in the diffusive limit: channels with low aspect ratios and small pores are subject to stronger shear (including recirculation bubbles for sufficiently large wall deformations), which acts to further enhance long-time dispersion.

Many biological transport processes have Péclet numbers in the range of Pe ∼ 10⁻²–10², and at the heart of all of these transport phenomena is the competition between entropic trapping and flow induced enhanced diffusion that is captured by our effective model. For example, the unicellular organism P. polycephalum has a complex network of veins that are used to deliver nutrients. It achieves this by generating contracting peristaltic waves that result in periodic shuttle flows. It has been reported that these flows can result in a two orders of magnitude increase in the effective diffusion coefficient compared to the molecular diffusivity of nutrient molecules.^51,52 Another relevant example is peristaltic transport in the small intestine, which occurs in the advection dominated limit with Pe ∼ 10⁶. It has been shown that the flows generated due to the rhythmic contraction of the walls in this system are crucial in determining the population density and transport of microbes that play important roles in the health of the host.²⁷ Consistent with the idea of Taylor dispersion, typical measurements suggest that these flows can enhance the bare molecular diffusivity by three orders of magnitude.

It is also interesting to note that there is a simple connection between the peristaltic flow and the classic problem of “Taylor’s swimming sheet.”^53–55 Taylor’s swimming sheet attains a non-zero swimming speed using traveling waves on its surface in contrast to peristalsis where the walls are not allowed to drift. There has been a considerable interest in the study of effective tracer diffusion and advective mixing in dilute bacterial suspensions.^56,57 We have carried out leading-order asymptotics to study dispersion by a swimming sheet in the vicinity of a wall, and the results are similar to those presented in this paper. The problem can be entitled in good humor as “Taylor dispersion by Taylor’s swimming sheet.”

Several other problems related to transport and dispersion in peristaltic flow remain to be explored. Many biological fluids involved in peristalsis are non-Newtonian in nature.⁴² Under long wavelength approximations, the leading order velocity field of a second-order fluid is identical to that of a Newtonian fluid,⁵⁸ and the asymptotic dispersivity is, therefore, also unchanged. However, the effects of viscoelasticity could be considered by including higher order corrections to the velocity field or employing a numerical solution for the velocity and B fields. Future studies may also address the geometric shape optimization of peristaltic waveforms for maximizing dispersion in both Newtonian and non-Newtonian fluids.^59,60 Finally, it is worth pointing out that the extension of this model for relevant biomedical applications comes with new challenges. For example, it is known that the enhancement in diffusion of passive tracers in blood capillaries is not only driven by macroscopic flows but is also significantly altered by micrometer-scale flow disturbances generated by individual red blood cells. Developing models that capture these subtleties is an area of active interest.⁶¹ 

The data that support the findings of this study are available within the article.

D.S. gratefully acknowledges funding from the National Science Foundation Grant No. CBET–1934199.

We validate the theory by carrying out Brownian dynamics simulations. The dimensionless Langevin equation corresponding to Eq. (18) in the moving frame is discretized as

where Δt is the time step and wⁿ(t) is a Gaussian random vector with zero mean and unit variance. The no-flux condition at the walls is imposed through the method of specular reflections.⁶² We typically carry out the simulations over an ensemble of N = 10⁵ particles with a time step of Δt = 10⁻⁶. The particles are released at R⁰ = 0, and we collect statistics at long times. The first two moments of the particle position are

The effective dispersivity then follows as