In biological transport mechanisms such as insect respiration and renal filtration, fluid travels along a leaky channel allowing material exchange with systems exterior to the channel. The channels in these systems may undergo peristaltic pumping which is thought to enhance the material exchange. To date, little analytic work has been done to study the effect of pumping on material extraction across the channel walls. In this paper, we examine a fluid extraction model in which fluid flowing through a leaky channel is exchanged with fluid in a reservoir. The channel walls are allowed to contract and expand uniformly, simulating a pumping mechanism. In order to efficiently determine solutions of the model, we derive a formal power series solution for the Stokes equations in a finite channel with uniformly contracting/expanding permeable walls. This flow has been well studied in the case in which the normal velocity at the channel walls is proportional to the wall velocity. In contrast we do not assume flow that is proportional to the wall velocity, but flow that is driven by hydrostatic pressure, and we use Darcy’s law to close our system for normal wall velocity. We incorporate our flow solution into a model that tracks the material pressure exterior to the channel. We use this model to examine flux across the channel-reservoir barrier and demonstrate that pumping can either enhance or impede fluid extraction across channel walls. We find that associated with each set of physical flow and pumping parameters, there are optimal reservoir conditions that maximize the amount of material flowing from the channel into the reservoir.

## I. INTRODUCTION

Flow through permeable channels at small Reynolds number occurs in a variety of physical and biological applications. Examples include air flow along trachea in an insect’s body^{1,2} or urine flow along collecting ducts in a mammal’s kidney.^{3} Flow may be driven by a variety of mechanisms; one well known mechanism is via muscle contractions that can dynamically deform the cross sectional area of transporting channels. In insect respiration, it is experimentally understood that pumping of the trachea, driven by thoracic or dorsal ventral muscle contractions, enhances oxygen intake and carbon dioxide output.^{2} However, it is unclear if the mechanism behind the enhanced respiration arises primarily from the increased advection, or from secondary pressure-driven flows across the tracheal membranes. In renal filtration a similar question persists: it remains controversial to what degree pumping, caused by muscle contractions along the pelvic wall, augments the molarity of waste products and salts in urine.^{3} The idea that pumping can lead to enhanced advection of transported fluids has been well studied.^{4–6} In the case of permeable channels, it is less well understood how pumping can alter material exchange across the channel walls. Such an understanding is crucial if one is to understand the mechanisms underlying enhanced respiration in the insect or filtration in the kidney.

To study this problem, we consider a simple filtration model in which fluid is exchanged between a channel and a reservoir (or interstitium). The problem of fluid flow through permeable channels has been studied extensively over the past 60 years. A common physical assumption is that the region exterior to the channel has a constant and fixed pressure, and fluid flow across the channel wall is driven by the hydrostatic pressure difference, typically modeled linearly with Darcy’s law. The seminal work on this problem was performed by Regirer,^{7} who determined an asymptotic solution for flow through a channel with low permeability and low Reynolds number. Subsequently, the problem of flow across channels driven by Darcy’s law has been expanded in a number of ways.^{8–15} A related problem considering prescribed cross-channel velocities in pipes and channels has also been the subject of significant investigation, beginning in 1953 with Berman,^{16} and this problem has also been extended in a variety of ways.^{17–19} Whereas the above works assume static channels, other investigations have examined channels and cylinders with deforming walls as a function of time.^{20–28} All of the studies involving moving channel walls rely on suction/injection boundary conditions, where the flow velocity at the wall is assumed to be proportional to wall velocity and independent of the hydrostatic pressure. Despite the physical reality that flow across channels is driven by pressure, we know of no analytic work that examines walls in motion while considering pressure driven flow across the channel walls and avoids the assumption of suction/injection boundary conditions.

In the present work we determine an exact solution for Stokes flow through a channel with permeable and uniformly moving walls. To derive this solution, we utilize the memoryless nature of the Stokes equations. A wall in motion at any given point in time will be equivalent to prescribing a velocity at the wall coupled with the hydrostatically driven velocity field for a fixed channel width. We will solve this static problem and use it to construct stream lines for a wall in motion at any given point in time.

With a channel flow solution, we may then calculate the flux between the channel and reservoir. We make the assumption that fluid in the reservoir may also be exchanged through pathways other than the channel-reservoir barrier. In the current model, we assume that this second exchange of reservoir fluid happens at a rate proportional to the amount of fluid in the reservoir. Physically, such a model suggests either a leak in the reservoir leading outside of the considered system or some process that consumes the fluid such as evaporation, combustion, transport, or metabolism. To close the relationship between reservoir pressure, and the amount of material within the reservoir, we consider a model in which the reservoir pressure is proportional to the amount of fluid within the reservoir; we note that more complicated closures may be incorporated; however, do not consider these in the present work.

With these assumptions, we model the homogeneous reservoir pressure *p _{res}* as

where *γ* is the rate of extraction or metabolism, *j* is the material flux entering or leaving the reservoir across the channel, and *α* is a constant of proportionality relating the material flux to the rate of change of pressure within the reservoir. The parameter *α* may be thought of as a stiffness coefficient of the reservoir; one simple physical interpretation of *α* is the pressurization of a piston with a spring in an attempt to return the system to an equilibrium volume. As material enters the reservoir, the spring is compressed leading to an increased pressure. The material flux may be derived by solving a system of fluid flow equations within the channel and then calculating the fluid flux across the boundary. We assume that the channel has a uniform width that changes in time. Along with pressure differences between the reservoir and the inlet/outlet conditions on the channel, dynamic changes in the channel width will produce flux across the channel-reservoir barrier.

Having found a solution to the fluid flow equations, we then determine the flux between the reservoir and channel at any point in time and numerically solve our model for the reservoir pressure. We find that changes in the amplitude and frequency of pumping can have significant effects on the amount of extracted fluid; the effects are shown to depend upon the channel length, the properties of the reservoir, and the channel permeability. Furthermore, we determine that for any fixed set of parameters not including *α*, there is a unique value of *α* that optimizes material extraction; such a result may have impact both on understanding the biological systems mentioned above as well as on filtration applications.

## II. PROBLEM STATEMENT AND SOLUTION

We begin this work by determining a formal power series solution for the flow described in the introduction; below we will utilize this solution in our model for *p _{res}* (see Eq. (1)). In two dimensions, the problem of Stokes flow with permeable boundary conditions and moving walls can be stated as

where *p* is the pressure field; *u*, *v* are the axial and cross sectional velocity fields, respectively; *P*_{±L} are the pressure boundary conditions at *x* = ± *L*, respectively, with *y* = ± *r*(*t*); *r*(*t*) is the half-channel width as a function of time (uniform in space); and *p _{res}* is the (uniform in space) pressure in the region of space exterior to the channel, considered to be some reservoir of fluid (see Fig. 1). We let

*x*and

*y*denote the axial and cross sectional coordinates, respectively. The choice of Darcy’s law along with no slip boundary conditions at the channel walls has been previously justified in the case of isotropic boundaries with small permeability

^{13}and can be further justified in biological applications where water moves across epithelial cell membranes through water channels in the normal, but not tangential, direction resulting in an anisotropic permeability for the membrane wall.

*p*may also be set to zero without loss of generality. We note that we may change the boundary conditions from pressure to averaged flow conditions, for example,

_{res} Finally, we will assume that *r* is a periodic function in time with period *T*.

To compute the solution, we divide the problem into two simplified problems via the linearity of the Stokes equations. We separate the pressure-driven flow from the wall-driven flow by noting that the solution to Equations (2)-(5) may be written as a superposition of solutions to the system given the pressure driven flow,

and the wall driven flow,

where

The pressure driven flow (Equations (7)-(10)) has been well studied; an exact solution^{15} and many approximations for small permeability at non-zero Reynolds number can be found in the literature (see, for example, Refs. 11 and 13).

We will therefore restrict our focus to Equations (11)-(14). We begin by nondimensionalizing this system at a fixed time with nondimensional variables,

Dropping the tildes and the subscript specifying wall driven flow, we arrive at the nondimensional equations

where

As previously noted, the memoryless nature of Stokes flows allows us to simply determine a solution for any fixed value of time.

### A. Solving the nondimensionalized system

To solve Equations (19)-(22), we will develop a formal expansion, letting

in which the normal velocity condition for the *n*th term satisfies the Darcy relationship for the residual pressure of the (*n* − 1)th term. To start the recursive procedure, we let the *n* = 1 term satisfy the boundary condition given by the wall motion. In other words, the boundary conditions for the normal velocity in each term are

which are used as boundary conditions for the set of equations given by the system

for all *n* ∈ ℕ.

Below we will demonstrate that this recursive formula leads to a formal power series solution for (*p*, *u*, *v*). We will then investigate the effects of the parameters $A$ and ℓ have on the convergence properties of this expansion.

#### 1. Inductive procedure to determine the formal series expansion

The first order term (at *n* = 1) is equivalent to the zeroth order expansion of the Berman problem.^{16} It has the solution

The normal velocity boundary condition for the second element of the expansion is then

We note that we have a constant plus a quadratic term in *x*. Due to the linearity of the Stokes equations, we can separate the two terms on the right hand side of Equation (32) and note that the constant boundary velocity term, $3A\u21132/2$, already gives rise to a solution promotional to the Berman solution presented in Equations (29)-(31). As shown below, the *O*(*x*^{2}) term found in Equation (32) will result in a residual pressure that is a fourth order polynomial in *x* with exclusively even coefficients. More generally, we will show that an *n*th degree polynomial of *x*^{2} for the *v _{n}* velocity at the boundary will lead to a residual pressure that is an (

*n*+ 1)th degree polynomial of

*x*

^{2}.

*If Equations (26)-(28) are closed with a normal wall velocity prescribed by an*

*n*

*th-degree polynomial of*

*x*

^{2}

*,*

*then the pressure*

*p*

_{n}*at the boundary will be an*(

*n*+ 1)

*th degree polynomial of*

*x*

^{2}

*,*

We present a proof of the proposition in Appendix A. Within the proof is an algorithm for constructing the higher order system. Because the pressure *p _{n}* is the boundary condition for the normal velocity

*v*

_{n+1}, this proposition allows us to recursively construct a formal power series solution for Equations (19)-(22), provided that the expansion converges. In Subsection II A 2, we will examine the convergence properties of this expansion conditioned on the nondimensional parameters $A$ and ℓ.

#### 2. Convergence

To demonstrate the convergence of the series expansion, we seek an upper bound on the radius of convergence. First we note that the potential function for the flow at the *n*th expansion term, *ψ _{n}*, may be written in terms of the highest order term in

*x*, given by Equation (A13) as

Letting the potential function of the solution be $\Psi =\u2211n=1\u221e\psi n$, we use the ratio test to find an upper bound on the radius of convergence in *x*, and find that $|x|<A\u22121/2$, which provides an upper bound for ℓ, the non-dimensional axial domain length, based on $A$.

Next, we analyze the convergence numerically by examining a variety of values for $A$ and ℓ. Due to the criteria of the ratio test that $A\u21132<1$, we plot the convergence regimes parametrized by $log(A)$ and $A\u21132$. There will be no error between the approximation $\Psi N=\u2211i=1N\psi n$ and the exact solution for all equations except the normal velocity boundary condition, presented in Equation (22). Thus we define the error to be ‖*V _{N}* − (1 +

*P*)‖

_{N}_{∞}, where $VN=\u2211n=1Nvn$ and $PN=\u2211n=1Npn$. We define the relative error to be the error normalized by the sup norm of the average between

*P*+ 1 and

*V*given as

We vary $A$ between 10^{−5} and 1 on a log scale and $A\u21132$ between 0.01 and 1. We increase *N* until either the relative error is below a threshold (set to be 10^{−4}) or until *N* = 25. We record the regions of convergence along with the depth of recursion needed to find convergence in Figure 2(a) and display the rate of convergence in Figure 2(b). We note that for Figure 2(b), we do not consider whether or not the error threshold has been achieved but instead report the order of the convergence.

We find that the upper bound of $A\u21132<1$ is not sharp over all parameter regimes. Indeed, the error is increasing between truncation value *N* = 24 and *N* = 25 for all values of $A\u21132>0.56$. In the current parameter runs, we have found an approximate upper bound on $A$ of 10^{−5/8}. We did not numerically investigate a sharper restriction as permeabilities greater than $A=10\u22125/8$ require an exceedingly short channel, with channel length less than one fifth of the channel width, in order for the recursive solution to converge. Thus the resulting formal solution is unlikely to be of any aid in modeling physical systems in this parameter regime. Our numerical results also reveal that the recursive technique will not converge when the length is too long or the channel width is too narrow by bounding ℓ for fixed $A$.

### B. Solution properties

The solution demonstrates that the assumption of suction/injection boundary conditions breaks down as the length of the channel grows; see Appendix B. This effect is independent of the non-dimensional permeability $A$. We also note that much of the previous work has been done in the context of examining self-similar solutions. There is nothing about our method of solution, nor our formal power series expansion, that suggests self-similarity, and yet we find that the velocity profiles are remarkably close to being self-similar and may indeed be so. This result is also presented in Appendix B. Despite the potential for a self-similar solution, we have not been able to use these assumptions and methods to construct a solution for Equations (19)-(22).

### C. Range of validity

We now consider the types of systems for which the above solution may be valid. We begin by remarking that as the channel becomes thinner, ℓ becomes larger for fixed length, and thus our solution will only be valid so long as $inft\u2208Rr(t)\u22650.56/\kappa \mu L3$. For fixed *κ* and *μ*, this sets an upper bound on the allowable length scales of the channel. In Ref. 13, the authors perform a brief literature review and find that $A$ is typically in the range of 10^{−6}–10^{−14} nondimensional permeability. This means that the maximal channel to tube length ratio ($maxt\u2113=maxtL/r(t)$) must be somewhere in the range of 7.5 × 10^{2}–7.5 × 10^{6} depending on the value of $A$. This implies that the solution will break down for fully compressed channels for which *r*(*t*) → 0, and that there is a minimal width of compression, below which the formal power series will no longer converge.

## III. AN AVERAGED FLOW MODEL FOR MATERIAL EXTRACTION

Having determined a solution for a specified range of parameters above, we now return to the original question: to what extent can pumping can alter averaged flow rates across the channel wall? We adopt a periodic, symmetric-in-time pumping strategy such that

In order to have a fair comparison between the pumping and non-pumping fluxes, we have chosen a pumping strategy that has a mean half-width $\u3008r(t)\u3009=r\u0304$. This will be compared to the case without pumping when the half-width is fixed at $r\u0304$. The flux in Equation (1) may now be described as

where *v* is the normal velocity from Equations (2)-(5). We subtract the speed of the wall in the interior of the integral of Equation (36) in order to obtain the flux across the wall, and then integrate to determine the total flow into the reservoir.

We will denote the pressure driven flow (the solution of Equations (7)-(10)) as *v _{press}*(

*x*) =

*v*(

_{p}*x*,

*r*(

*t*)) and the wall motion-driven flow (the solution of Equations (11)-(14)) as

*v*(

_{pump}*x*) =

*v*(

*x*,

*r*(

*t*)). From the solution to Equations (7)-(10),

^{15}we determine that the normal velocity at

*y*=

*r*(

*t*) from the pressure driven flow is given by

where *λ*(*t*) is the solution to the transcendental equation

We note that the leading order term in the results found in Refs. 7 and 13 are asymptotic approximations for Equation (37) in the limit of $A\u21920$, and the only difference in Equation (37) and the asymptotic approximations is that one sets $\lambda (t)\u22483A(t)$. Either the exact solution found in Ref. 15 or the asymptotic theory resulting from Refs. 7 and 13 may be used in the following exposition; we have tested both theories and found no difference in the results.

To nondimensionalize the system, we set

and find that the non-dimensional pressure driven flux to be

where

and the nondimensional constants are defined as

and the time varying nondimensional parameters are defined as

Having found the flux across the channel and reservoir interface for pressure driven flows, we are now ready to determine the averaged flux for our given pumping strategy. Before turning to study the model for *p _{res}* presented in Equation (1), we will begin by examining an intermediate result in which

*p*is a fixed value,

_{res}*p*(

_{res}*t*) =

*P*.

_{res}### A. Intermediate model: *p*_{res}(*t*) = *P*_{res}

_{res}

_{res}

In the case when *p _{res}*(

*t*) =

*P*, we seek to determine the relative change in the average flux when we add pumping into the system. We are thus interested in the quantity,

_{res} where *j*_{0} is the flux in the absence of pumping and *j* is the flux with pumping. One may readily see that the averaged wall driven flux is zero. This may be seen by noting that the non-dimensional solution for the velocity at the wall velocity due to pumping is independent of the sign of the wall velocity; physically, this is a consequence of reversibility. The pressure driven flux, one the other hand, is not symmetric as can be seen from Equation (40); physically, this corresponds to the fact that halving a channel width restricts the fluid flux more significantly than doubling the channel width enhances the fluid flux.

Because Equation (44) acts as a linear operator on the flux, we ignore pumping and simply consider the pressure driven flux. Plugging Equation (40) into Equation (44) yields

where $j\u0303press0$ is the pressure driven flux in the absence of pumping. We have not found a closed form expression of this integral and therefore proceed numerically. We consider the restriction in parameter space by bounding 0 < Θ < 0.5. As mentioned above, we consider realistic values for the nondimensional permeability, confining $A\u0304\u2208[10\u221214,10\u22126]$. These two constraints determine a maximal value for $\u2113\u0304$ of 265, to ensure we remain within the radius of convergence of the formal power series for *v*_{pump}. With these constraints in mind, we proceed by setting $A\u0304=10\u22126$, $\u2113\u0304=200$, and vary |Θ| from 0.1 to 0.5. Our result will be independent of the choice of Ω. We numerically evaluate Equation (45) and present the integrand along with the averaged flux as a function of Ω in Figures 3(a) and 3(b), respectively. We find that in fixing the reservoir pressure, pumping is exclusively able to impede the flow independent of its direction. Next we will investigate the more complicate reservoir dynamics found in Equation (1) and demonstrate that such dynamics may actually lead to enhanced flows; in fact in this new model we will demonstrate that even when there is no directional flow in the absence of pumping, pumping can generate directional flow.

### B. A model with dynamic reservoir pressure

Employing the same nondimensionalization as found in Equation (39), we nondimensionalize Equation (1) and find

where *v _{pump}* is the (nondimensional) solution for Equations (19)-(22) parametrized by $A(t\u0303)$ and $\u2113(t\u0303)$. The nondimensional constants are the same as above (see Equations (42) and (43)).

We thus arrive at a linear first order differential equation characterized by six nondimensional parameters (Equation (42)). We may express the exact solution of this equation in terms of an integral equation

where

By Floquet theory, a time-periodic solution exists, and since the coefficients are smooth, this solution is unique up to transient behavior. Below we examine the effect of selected parameters on the mean value of the reservoir pressure and we show that pumping can indeed enhance the average material extraction.

#### 1. Generating directional flow

We begin by considering the case in which there is no flow in the absence of pumping. Mathematically, this corresponds to setting Ξ = 0. In this reduced space, we consider two thought experiments.

We first consider the limit of small Υ. Physically, this corresponds to a compliant reservoir in which incoming fluid has little effect on the reservoir pressure. We note that the coefficients on the differential equation have period 2*π*/Ω. Suppose then that we have a periodic solution and expand the solution in the Fourier basis, setting Υ = 0. The zeroth mode gives the mean value of the solution which is only influenced by the zeroth mode of $f(t\u0303)$. We note that average value of $f(t\u0303)$ is the amount of fluid that enters and exits the reservoir after a contraction and expansion due solely to the motion of the wall and therefore has mean zero due to the reversibility of the Stokes equations. Thus periodic solutions for the reservoir pressure have zero mean. Physically, such a model supposes that pressure driven flows are dominated by motion driven and metabolic flows. Due to the reversibility of the Stokes equations, this leads to zero net flow. This explains the case in which reservoir does not change its pressure with respect to volume.

Next we consider the limit of large Υ. Physically, this corresponds to a stiff reservoir in which incoming fluid must either immediately leave through the metabolic process or be reabsorbed into the channel. In this case, $f(t\u0303)$ is small compared to $g(t\u0303)\u22121$, i.e., ΘΩ/Υ = *ϵ* ≪ 1 (since Ξ = 0). We have determined in Sec. II that $f(t\u0303)$ is absolutely bounded within its radius of convergence, and further note that $g(t\u0303)>1$ for all time, which implies that the first term of the convolution is bounded by 1. Letting *M* be the absolute bound for $f(t\u0303)$, we conclude that the convolution found in Equation (47) is absolutely bounded by *ϵM*. Furthermore, the first term of Equation (47) decays to zero in long time. It follows that the reservoir pressure goes to zero as *ϵ* → 0 and $t\u0303\u2192\u221e$. Physically, this corresponds to the situation in which pressure driven flow dominates the dynamics. In this case, the system will not exchange material due to the growing energetic cost of pressurizing the reservoir.

If any flow direction can be produced, it must be due to a balance between pressure driven and motion driven fluxes, which corresponds to intermediate values of Υ with the other parameters fixed. To numerically investigate the effect of intermediate values of Υ, we consider the restriction in parameter space of bounding 0 ≤ Θ ≤ 0.5. We proceed by setting Ω = 10, $A\u0304=10\u22126$, $\u2113\u0304=200$, and vary |Θ| from 0.1 to 0.5, and Υ from 0 to 0.05. We then determine the long time average of $p\u0303res(t\u0303)$,

The reason we seek to determine the long time average of $p\u0303res(t\u0303)$ is that it is the average metabolic material flux. Additionally it is proportional to the pressure driven flux. Thus if $\u3008p\u0303res(t\u0303)\u3009$ is positive, we have directional flow from the channel into the metabolic process, and if it is negative, we have flow into the channel from the reservoir reversing the metabolic pathway. We solve the system using Mathematica’s built-in NDSolve routine.^{29} For the remainder of the paper, we truncate the formal series solution for *v _{pump}* at

*N*= 5 since (i) the relative two norm is below 0.5% and (ii) we have noticed no difference in the results presented below in adding terms to the series expansion. Results are shown in Figure 10. For each value of Θ, an optimal value of Υ for material extraction can be found between 0.016 and 0.021, which implies parameter regimes can be identified in which pumping generates maximal enhanced material extraction (Fig. 4). In other words, there is an Υ that maximizes material extraction. We further note that the location of this maximal value is fairly insensitive to Θ.

We continue a sensitivity analysis in Appendix C 1 and find that Υ is also insensitive to changes in $\u2113\u0304$ and $A\u0304$, and that it grows with Ω. The amount of extracted material grows with Ω, but is bounded from above. Thus, we conjecture that with Ξ = 0, no matter how we set the parameters $\u2113\u0304$, $A\u0304$, Ω, and Θ, we will find value for Υ that maximizes the amount material being extracted from the reservoir. Furthermore, we conclude that pumping is able to generate an average flow direction across the channel-reservoir interface.

Although we have numerically demonstrated the existence of directional flows, we have yet to explain the mechanism driving this flow. In Appendix C 2, we hypothesize several potential mechanisms and determine that it is the phase difference between the pressure flux and the pumping flux that generates the directional flow. We also offer an explanation for why this phase communication occurs and note that the consequence for the physical model is that net flow will only be in the channel to reservoir direction (see Appendix C 2).

#### 2. The effect of Ξ

We conclude our parameter study by changing the average amount of material being extracted in the absence of pumping and determine the percent gain in material extraction when pumping is introduced into the system. Parametrically, this means that we change Ξ and note that under no pumping (either Ω = 0 or Θ = 0), $p\u0303res$ has a stable steady state value given by

Similar to the above study in Section III B 1, we set $A\u0304=10\u22126$, $\u2113\u0304=200$, Ω = 10, Θ = 0.5 and vary Υ; for each value of Υ, we set Ξ so that $p\u0303resss=1$. Simulation results demonstrate that for smaller values of Υ, material extraction is impeded by roughly 4% (see Figure 5). To identify the mechanism of this impedance, we first note that for $\u2113\u0304$ and $A\u0304$ fixed, the average pressure is expected to converge to $p\u0303resss=1$ in the limit that either Θ or Ω approach zero. We numerically test this idea by decreasing both Ω and Θ and present our results in Table I. We find that the reservoir pressure does *not* approach the steady state value in the limit of Ω → 0 and instead decreases until it converges to some value less than 1. As Θ → 0 however, we do see $p\u0303resss\u21921$.

Ω | 10 | 1 | 0.1 | 0.01 |

$\u3008p\u0303res\u3009$ | 0.9631 | 0.9623 | 0.9614 | 0.9614 |

Θ | 0.4 | 0.3 | 0.2 | 0.1 |

$\u3008p\u0303res\u3009$ | 0.9791 | 0.9894 | 0.9956 | 0.9989 |

Ω | 10 | 1 | 0.1 | 0.01 |

$\u3008p\u0303res\u3009$ | 0.9631 | 0.9623 | 0.9614 | 0.9614 |

Θ | 0.4 | 0.3 | 0.2 | 0.1 |

$\u3008p\u0303res\u3009$ | 0.9791 | 0.9894 | 0.9956 | 0.9989 |

In the case that Ω → 0, the contribution of the mass flux due to pumping becomes small when compared to the mass flux driven by the hydrostatic pressure. This suggests that the third term on the righthand side of Equation (46) becomes small. Assuming that the contribution of flux from pumping remains bounded (as it will within the radius of convergence), we can write the solution to Equation (46) with Ξ ≠ 0 as

where

To study this solution, we first ignore the transient term by setting $p\u0303res(0)=0$ and neglect the small *O*(Ω). Next we set Υ = 0 to ensure that we examine the regime of parameter space in which we expect material extraction to be impeded by pumping (see Figure 5). Finally, since we will examine the small limit of Ω, we set $t\u0303=2\pi \tau /\Omega $ where *τ* ∈ [1, 2] so that we may examine flow over a full period of pumping. We do not choose *τ* ∈ [0, 1] in order to avoid considering the transient behavior for small values of $t\u0303$. We find that as Ω → 0,

Figure 6 shows $p\u0303res$ as a function of *τ* for several values of Θ, $A\u0304$, and $\u2113\u0304$. We find that the reservoir pressure decreases when the channel is compressed and vice versa; however, this effect is asymmetric. Indeed, this asymmetry enables the overall average to decrease in the mean pressure reservoir for low values of Ω. As demonstrated in Figure 5, this effect is seen even when Ω = 10. We note that this is the identical effect found in the intermediate model in which we fixed the reservoir pressure.

## IV. DISCUSSION

We have developed a model for material extraction and demonstrated both qualitatively and quantitatively the mechanism by which uniform pumping enhances material extraction from a channel with permeable walls. Furthermore, we have provided evidence of optimal conditions for material extraction in terms of the ratio between reservoir stiffness and material depletion. To determine the flux across the channel-reservoir barrier, we have derived a formal series expansion for Stokes flow in a channel with permeable and uniformly moving walls. Although there has been a great deal of work examining this problem (see, for example, Refs. 20–28), we have examined a solution in which pressure gradient drives flow across the boundary rather than suction/injection boundary conditions. To our knowledge, our work is the first to examine this boundary condition in the context of expanding and contracting walls. We have investigated the convergence regimes of this series expansion and found analytic bounds for the parameter domain of convergence, and have numerically determined the convergence regimes with respect to the nondimensional permeability and channel length.

The mechanism that facilitates material extraction is fairly deep and is explored in Appendix C 2. The simple explanation for this occurrence is that when the pressure in the reservoir is the highest, the pressure driven flow across the reservoir channel interface is slow because the channel width is reduced. Mathematically, we find that it is the natural phase difference between the reservoir pressure along with the fluxes due to pressure and wall motion that facilitate the existence of flow; we note however that one has no control over the phase difference between the coefficients on the model presented in Equation (1). Because of this lack of control, we have found that pumping will only generate averaged flow leading from the channel into the reservoir. It would be interesting to determine if there exists a pumping strategy (*r*(*t*)) that would lead to flow entering the channel in a future study.

The simplicity and limitations of the present model when applied to biological processes such as insect respiration or renal filtration must be acknowledged. One discrepancy is that while the present model represents uniformly moving walls, compression of the biological channel typically occurs at a discrete location (insect respiration^{4}) or in a wave like pattern (renal filtration^{3}). Second, exchange across the channel-reservoir barrier is typically not as simple as a simple homogeneous fluid. In insect respiration, air may be thought of as a mixture of oxygen and unusable elements such as nitrogen and carbon dioxide. To model this distinction, a phase field model may be more suitable for the air flow, and in addition only oxygen should be depleted (and converted into waste) within the reservoir. However, if the permeabilities for the waste and oxygen are equivalent, then the phase field model would decouple from the flow equations, and the work presented here would be a first step in solving such a model. In renal filtration, not only water is being filtered, but ions as well. Thus osmotic pressures must be taken into account when considering fluid exchange across the channel-reservoir barrier. Despite these limitations, having obtained a promising first result on the effect of pumping on fluid extraction, we plan to investigate these more realistic models in future work.

## Acknowledgments

This research was supported by the National Science Foundation (NSF) Grant No. DMS1263995, to A. Layton; Research Training Groups Grant No. DMS0943760 to the Mathematics Department at Duke University; KI-Net NSF Research Network in the Mathematical Sciences Grant No. 1107291, and NSF Grant No. DMS1514826 to Jian-Guo Liu.

### APPENDIX A: PROOF OF PROPOSITION 1

We will prove Proposition 1 by by induction. In the proof, we present a simple algorithm for constructing the higher order systems.

*If Equations (26)-(28) are closed with a normal wall velocity prescribed by an*

*nth-degree polynomial of*

*x*

^{2}

*,*

*then the pressure*

*p*

_{n}*at the boundary will be an*(

*n*+ 1)

*th degree polynomial of*

*x*

^{2}

*,*

To prove the proposition by induction, we first note that the base case has already been satisfied (Equations (29)-(31)). Supposing the inductive hypothesis holds true, we may note that each individual term in the polynomial for *v _{n}* with

*j*<

*n*− 1 leads to a residual pressure that is a polynomial of degree at most (

*n*− 1) of

*x*

^{2}at the boundary. Therefore we must demonstrate that

leads to a residual pressure that is a polynomial of degree at most *n* of *x*^{2} at the boundary. We set *a _{n}* = 1 without loss of generality, due to the linearity of the equations.

Let *ψ _{n}* be the flow potential such that

*u*= ∂

_{n}_{y}

*ψ*and

_{n}*v*= − ∂

_{n}_{x}

*ψ*. Due to the symmetry of the boundary conditions on the pressure and velocity at the boundaries, we assume that the flow profile

_{n}*u*is odd in

_{n}*x*and even in

*y*, that

*v*is even in

_{n}*x*and odd in

*y*, and that

*p*is even in both

_{n}*x*and

*y*. Next, let us suppose there exists a finite power series such that

This ansatz is natural in that it (i) ensures the symmetry conditions listed above and that it (ii) may satisfy the boundary conditions. We note that should this ansatz give rise to a solution, then the pressure will be given by the indefinite integral

with an unknown constant of integration *f* that is a function of *y* alone. By directly integrating Equation (A2) and then setting *y* = 1, one obtains that the pressure at the boundary is a polynomial in *x* of degree 2(*n* + 1). The symmetric boundary condition at *y* = − 1 is satisfied due to the symmetry of the ansatz presented in Equation (A2).

We then check condition (ii), which is equivalent to showing that the ansatz leads to a solution to the system given in Equations (26)-(28), or that

The only unknowns left in the system are the *c _{ij}*’s of which there are (

*n*+ 3)(

*n*+ 2)/2. Plugging

*ψ*into the biharmonic equation and equating powers

_{n}*x*and

*y*results in (

*n*+ 1)

*n*/2 equations described by

for *j* ∈ {1, *n*} and *i* ∈ {1, *n* + 1 − *j*} due to the fact that when *i* = 0 or *j* = 0, the linear terms in *x* and *y* (respectively) vanish upon the first application of the Laplacian. We also note that $Pkn$ represents the *k*-permutations of *n*. The normal velocity condition and no slip condition require that

for *i* ∈ {0, *n* + 1} and where *δ _{ij}* denotes the Kronecker-

*δ*function. This system leads to an additional 2(

*n*+ 2) equations. According to the above calculations, there appears to be one more equation than unknown, hinting that the system may be over determined. However, there is a trivial redundancy when

*i*=

*n*+ 1 in Equations (A8) and (A9) which give, respectively,

Thus we have reduced the system to a square matrix equation. Next we will show the nonsingularity of the matrix by developing a recursion formula to directly solve the system. To do this, we note from above that *c*_{(n+1)0} = 0 and that we have (*n* + 2 − *i*) unknowns with two equations for every fixed value of *i* using Equations (A8) and (A9). We have already determined the case where *i* = *n* + 1. For *i* = *n*, we obtain the 2 × 2 dimensional system

which has solution

Proceeding recursively, we assume that we have found *c _{ij}* for all values of

*i*′ greater than some fixed

*i*and

*j*= {1, …,

*n*+ 1 −

*i*′} for each value of

*i*′ >

*i*. We seek to demonstrate that Equation (A7) can be used to determine

*c*for

_{ij}*j*≥ 2, thereby reducing the (

*n*+ 2 −

*i*) to two parameters (namely,

*c*

_{i0}and

*c*

_{i1}). From here, we note that Equations (A8) and (A9) are reduced to

where

are known quantities. This system is linearly independent. Thus the entire system is invertible and a unique polynomial exists. The remaining claim is straightforward as we have Equation (A7) for all values of *i* ≥ 1, and thus we can find *c _{ij}* (for 1 <

*j*<

*n*+ 1 −

*i*) by setting it as a function of

*c*

_{(i+1)(j−1)}and

*c*

_{(i+2)(j−2)},

### APPENDIX B: SOLUTION PROPERTIES

Having determined a formal series solution and numerically investigated convergence properties in Section II A, we now investigate the properties of our solution in the context of previous work (see below). We begin by analyzing the assumption of suction/injection boundary conditions. We continue by determining the proximity of our solution to a self-similar solution in the axial direction, as much of the work in this field has investigated solutions with this property.

#### 1. A comparison with the assumption of suction/injection boundary conditions

As previously noted, there has been a considerable amount of work examining solutions with the assumption of suction/injection boundary conditions (see, for example, Refs. 20–28). These boundary conditions assume that the velocity at the wall is proportional, but not equal, to the wall velocity at all points along the tube. That is, the normal boundary velocity (Eq. (5)) is given by

where the permeance coefficient *c* represents the ratio between the wall speed and the flow velocity at the wall. In nondimensional units, this boundary condition becomes

It is natural to ask how our model, driven by a hydrostatic pressure difference across the channel walls, compares to a model assuming suction/injection boundary conditions. To analyze this, we ignore the pressure driven flow, use our recursive solution to Equations (19)-(22), and approximate

We then determine the relative variation away from the assumption of suction/injection boundary conditions as

*σ _{si}* is constructed so that it is invariant under linear spatial transformations. As shown in Figure 7, the assumption of suction/injection boundary conditions breaks down for all values of the nondimensional permeability, $A$, as the ratio between channel length and width, ℓ, increases. Additionally, for values of $A<0.1$, the variation

*σ*is roughly fixed when a value for $A\u21132$ is fixed. To see this concretely, we overlay plots of

*σ*as a function of $A\u21132$ for values of $A\u2208[10\u22126,3.16\xd710\u22123]$ in Figure 8.

_{si}#### 2. Flow profile variation along the axial direction

Perhaps the most common technique in developing analytic theory of flow through permeable channels and pipes is to assume a self-similar profile along the axial direction. The present series expansion is not self-similar (see Equation (33)). Nevertheless, one may ask how the profiles of the field variables change with respect to the axial position.

We proceed both analytically and numerically. For the analytic aspect, we first take a different nondimensionalization of Equations (2)-(5) by letting

Setting the boundary pressures to be zero, we arrive at a new system of nondimensional equations (dropping tildes)

This system may be solved using an algorithm identical to that presented in Section II A, with the exception of an extra-power of $A$ in the expansion for *p* which is accounted for in the new nondimensional units. Thus we may now view the formal series expansion as an asymptotic expansion about small values for $A$. This is a reasonable restriction given that the series only converges for sufficiently small values of $A$. We then note that we may approximate our solution with the first two terms in the asymptotic expansion, where the zeroth order term, $O(A0)$, corresponds to Berman’s solution. This expansion predicts that

Although these approximations are, strictly speaking, not self-similar profiles, the deviation is small. This can be seen by considering the variance in the profile

Because the integrands are polynomials, we can analytically determine *σ*_{p}(*u*) and *σ*_{p}(*v*) for any given set of parameters $A$ and ℓ. In Figure 9, we plot max(*σ*_{p}(*u*), *σ*_{p}(*v*)) and see that with the leading order asymptotic theory, the profile variation stays below 5 × 10^{−3}. We repeat the analysis numerically for profiles that have converged below the error tolerance. To save on the numerical expense of directly calculating *σ*_{p}(*w*), we approximate it by setting $\beta =maxy(|w(x1,y)/w(x2,y)|)$ and consider the limit as *x*_{1} → ℓ and *x*_{2} → 0. Our choice for *β* ensures that computed value will be an upper bound for *σ _{p}*; our choice for

*x*

_{1}and

*x*

_{2}is made because numerical test cases suggest this will yield the greatest difference in the profile shape. We find that the profiles vary less as more terms in the series expansion are added. All profile variations for the numerical results were bounded below 5 × 10

^{−7}, which is a surprising result considering that the error tolerance on the numerical scheme was set to 10

^{−4}. This may explain the fluctuations in Figure 9(b). There is nothing

*a priori*to suggest the solution should form a self-similar profile. Nevertheless, we have shown that the profile may be well approximated by a self-similar profile.

### APPENDIX C: PROPERTIES OF MATERIAL EXTRACTION

#### 1. A parameter study on the fluid extraction model

To assess the sensitivity on the other parameters for the model presented in Section III B, we consider $\u3008p\u0303res\u3009$ as a function of Υ for Ω ∈ [10, 100], $\u2113\u0304\u2208[100,200]$, and $A\u0304\u2208[10\u22128,10\u22126]$ (see Figures 10(b)-10(d), respectively); additionally, Θ = 0.5, Ω = 10, $A\u0304=10\u22126$, and $\u2113\u0304=200$. Our results indicate that the optimal Υ is insensitive with respect to Θ and $A\u0304$, changes with respect to $\u2113\u0304$, and is highly sensitive to changes in Ω. This suggests that one ought to design a reservoir with a particular frequency in mind and then vary Θ to change the amount of material extraction.

Qualitatively similar results were obtained for other parameters. For example, we present the results obtained for Θ = 0.5, Ω = 0.1, $A\u0304=10\u22127$, $\u2113\u0304=500$, and Ξ = 0 in Figure 11. We have not found any set of parameters for which pumping extracts material away from the reservoir when Ξ = 0.

##### a. The dependency of $\u3008p\u0303res\u3009$ on $A\u0304$, $\u2113\u0304$, and Ω

Having determined the existence of an optimal value for Υ for several parameter ranges, we next inquire how changes in the other parameters alter the average amount of material extracted. Material extraction as a function of $A\u0304$, $\u2113\u0304$, and Ω is shown in Figure 12. By default $A\u0304=10\u22126$, $\u2113\u0304=200$, Ω = 10, Υ = 0.02, Θ = 0.3, and Ξ = 0. Our results indicate that the mean reservoir pressure increases with Ω, $A\u0304$, and $\u2113\u0304$. For Ω, the gain appears to be bounded from above , implying that increasing the rate of pumping is accompanied by diminishing returns beyond a certain frequency. Qualitatively similar results are obtained for other points in parameter space (not shown).

#### 2. The mechanism of flow

We turn next to determine the mechanism behind the generated flow presented in Section III B 1. Because the mechanism is not mathematically obvious, we adopt a scientific approach and visualize the qualitative form of the coefficients on the differential equation in Figure 13, in the case that Ξ = 0, $A\u0304=10\u22126$, $\u2113\u0304=200$, Θ = 0.5, Υ = 0.01, and Ω = 10. We pick these parameters case since we have seen that they will give us averaged directional flow. We already know that dramatic changes in the relative scaling between these coefficients will cause the directional flow to vanish, since such a change is equivalent to changing Υ and we know that $\u3008p\u0303res\u3009\u21920$ when Υ → 0 or Υ → ∞. Thus, without changing scales, we propose three hypotheses as to why we observe enhanced material extraction:

Pumping generates material extraction because of the skewness of the motion driven flux ($f(t\u0303)$).

Pumping generates material extraction because of the irreversibility in the pressure driven coefficient ($g(t\u0303)\u22121$). This is the same mechanism as the intermediate modeling with the opposite effect, namely enhanced flow rather than impeded flow.

Pumping generates material extraction because of the phase communication between coefficients.

We test hypothesis 1 by an approximate system in which the magnitude (in the supremum norm) is identical but in which we have control over the skew of *f*. Because we are no longer presenting a physical model, we change the differential equation to read

where

where *c*_{1} is some constant,

and *s* ∈ (−∞, ∞) is related to a measure of the skewness. We restrict the range of *s* ∈ (−10, 10) as the extremes of this range that are large enough to generate a sawtooth wave (*s* positive and large) and a time reversed sawtooth wave (*s* negative with large magnitude). We set *c*_{1} = 50 and vary *s*. We note that the material extraction was independent of the sign of *s*. Results are shown in Figure 14. We find that we achieve material extraction independent of the skew, although the amount of material extraction does change. It is interesting to note that a non-skewed function enhances the material extraction; however, the value of the skew is not the cause of the material extraction. Thus, the first hypothesis has been eliminated as a mechanism of directional flow.

To test hypothesis 2, we construct a function in which we can control asymmetry of the *g* function about its mean. To do this, we consider a second simplified model

where

*c*_{2} is some constant,

and *r* ∈ (0, 1) is a parameter that controls the symmetry of *G* about its mean *c*_{2}. When *r* is small, we will recover a sine wave and when it is large, we will have a one sided cusp. We numerically solve this system for the parameters mentioned above and varying *r*. We fix *c*_{2} = 1. We find that the averaged flow profile is nearly independent when changing the symmetry about the mean of the *G* function, which falsifies the second hypothesis (Fig. 15).

In testing hypothesis 3, we again consider a simplified model given by

In such a model, we have course grained out many of the important physical features, but have retained a sense of relative importance between “pressure” driven (Υsin(Ω*t*)) and “motion” driven flows with the Υ parameter, along with the time ratio between material loss and pumping frequency with the Ω parameter. In this case Υ ∈ [0, 1] to ensure the “pumping” term is sign definite as it is in the physical model. We have also introduced a phase shift parameter *ϕ* which we will adjust to test the hypothesis. We begin by fixing Ω = {0.1, 1, 10} and Υ = 1. We find that by altering the phase difference, we are able to completely reverse the flow direction (see Figure 16(a)). There is thus strong evidence that the communication between phases of the pressure driven and velocity driven flow enables pumping to generate averaged directional flow. Finally, we note that the frequency may give us control over the degree of magnified flow; however, we do not observe flow reversal upon changing the pumping frequency.

To verify this, we return to the original nondimensionalized model (Equation (46)), but add an artificial phase shift to the *g* function presented in Equation (48). The test model becomes

where *ψ* ∈ [0, 2*π*/Ω). Similar to the test on the model in Equation (C7), we examine the model for a variety of pumping frequencies. We present the results in Figure 16(b) and find again that frequency changes the location of maximal and minimal flow extraction, and also that in the absence of a phase shift, changes in frequency are not able to reverse the flow direction.

To conclude, we seek to gain a stronger physical intuition for the mechanism behind the directional flow. We consider the simplified model from Equation (C7) with Υ = 1 and Ω is small so that the metabolic pathway dominates the recalibration of the reservoir pressure at any given time. Our physical interpretation is that the flux due to metabolism is −*y*, the flux due to interfacial motion is due to −(1 + Υsin(Ω*t* + *ψ*)) *y*, and the flux due to wall motion is due to the final term sin(Ω*t*). When the system is in phase (*ψ* = 0), we note that when the most fluid is being input into the system (*t* = (2*πn* + *π*/2)/Ω), the primary mechanism for repressurization is through back flow into the channel, and thus much of this fluid is not metabolized. When the most fluid is being taken out of the system (*t* = (2*πn* + 3*π*/2)/Ω), the dominant method of repressurization is through the reverse metabolic pathway. Because of this, we have a pump that reverses the metabolic pathway. When Ω grows the effect is the same, but we must consider transient effects of the reservoir and thus find a phase shift in the location of the maximum and minimum extraction points.