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.

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 body1,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 pres as

ṗres=γpres+αj(t,pres),
(1)

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.

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 pres (see Eq. (1)). In two dimensions, the problem of Stokes flow with permeable boundary conditions and moving walls can be stated as

p(L,±r(t))=PL,p(L,±r(t))=PL,
(2)
p=μΔu,u=0,
(3)
u(x,±r(t))=0,
(4)
v(x,±r(t))=±[κ(p(x,±r(t))pres)+ṙ(t)],
(5)

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

r(t)r(t)u(±L,y)dy=Ū±L.
(6)

Finally, we will assume that r is a periodic function in time with period T.

FIG. 1.

A schematic of the physical set up. The channel walls are moving toward each other in the direction of the arrows and the pressure is fixed at the corner positions. The pressure exterior to the channel is uniform.

FIG. 1.

A schematic of the physical set up. The channel walls are moving toward each other in the direction of the arrows and the pressure is fixed at the corner positions. The pressure exterior to the channel is uniform.

Close modal

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,

pp(L,±r(t))=PL,p(L,±r(t))=PL,
(7)
pp=μΔup,up=0,
(8)
up(x,±r(t))=0,
(9)
vp(x,±r(t))=±[κ(pp(x,±r(t))pres)],
(10)

and the wall driven flow,

pw(L,±r(t))=0,pw(L,±r(t))=0,
(11)
pw=μΔuw,uw=0,
(12)
uw(x,±r(t))=0,
(13)
vw(x,±r(t))=±[κ(pw(x,±r(t)))+ṙ(t)],
(14)

where

puv=ppupvp+pwuwvw.
(15)

The pressure driven flow (Equations (7)-(10)) has been well studied; an exact solution15 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,

x̃=x/r(t),ỹ=y/r(t),
(16)
ũ=u/ṙ(t),ṽ=v/ṙ(t),
(17)
p̃=pκ/ṙ(t).
(18)

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

p(,±1)=0,p(,±1)=0,
(19)
p=AΔu,u=0,
(20)
u(x,±1)=0,
(21)
v(x,±1)=±(p(x,±1)+1),
(22)

where

=Lr(t),A=μκr(t).
(23)

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

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

puv=n=1pnunvn
(24)

in which the normal velocity condition for the nth 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

vn(x,±1)={±1,n=1±pn1(x,±1),n>1,
(25)

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

pn(,±1)=0,pn(,±1)=0,
(26)
pn=AΔun,un=0,
(27)
un(x,±1)=0,
(28)

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

u1(x,y)=3x2y21,
(29)
v1(x,y)=12(3yy3),
(30)
p1(x,y)=3A2(x2y22+1).
(31)

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

v2(x,±1)=3A2(x22).
(32)

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, 3A2/2, already gives rise to a solution promotional to the Berman solution presented in Equations (29)-(31). As shown below, the O(x2) 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 nth degree polynomial of x2 for the vn velocity at the boundary will lead to a residual pressure that is an (n + 1)th degree polynomial of x2.

Proposition 1.
If Equations (26)-(28) are closed with a normal wall velocity prescribed by annth-degree polynomial ofx2,
vn(x,±1)=±j=0najx2j,
then the pressurepnat the boundary will be an (n + 1)th degree polynomial ofx2,
pn(x,±1)=±j=0n+1bjx2j.

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 pn is the boundary condition for the normal velocity vn+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 nth expansion term, ψn, may be written in terms of the highest order term in x, given by Equation (A13) as

ψn(x,y)=Any3x2n+13yx2n+12(2n+1)+o(x2n1).
(33)

Letting the potential function of the solution be Ψ=n=1ψn, we use the ratio test to find an upper bound on the radius of convergence in x, and find that |x|<A1/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 A2<1, we plot the convergence regimes parametrized by log(A) and A2. There will be no error between the approximation ΨN=i=1Nψ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 ‖VN − (1 + PN)‖, where VN=n=1Nvn and PN=n=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

ϵ(N)VN(1+PN)12(VN+(1+PN)).
(34)

We vary A between 10−5 and 1 on a log scale and A2 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.

FIG. 2.

The number of recursive steps needed to obtain a relative error of 10−4 is displayed (a). The maximum number of steps taken is limited to 25. We have only investigated parameters above the dashed line which displays A2=0.01. In order to determine the rate of convergence, and the existence of other parameters for which the recursive solution converges, we plot the rate of convergence given by a linear fit of log10(ϵ(N)) as a function of N. Regions colored white (and above the dashed line) signify diverging power series (b). We plot the line A2=1 which is an analytically derived upper bound of this product and A2=0.56 which is the numerically derived upper bound.

FIG. 2.

The number of recursive steps needed to obtain a relative error of 10−4 is displayed (a). The maximum number of steps taken is limited to 25. We have only investigated parameters above the dashed line which displays A2=0.01. In order to determine the rate of convergence, and the existence of other parameters for which the recursive solution converges, we plot the rate of convergence given by a linear fit of log10(ϵ(N)) as a function of N. Regions colored white (and above the dashed line) signify diverging power series (b). We plot the line A2=1 which is an analytically derived upper bound of this product and A2=0.56 which is the numerically derived upper bound.

Close modal

We find that the upper bound of A2<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 A2>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=105/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.

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

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 inftRr(t)0.56/κμ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=maxtL/r(t)) must be somewhere in the range of 7.5 × 102–7.5 × 106 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.

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

r(t)=r̄+Asin(ωt).
(35)

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 r(t)=r̄. This will be compared to the case without pumping when the half-width is fixed at r̄. The flux in Equation (1) may now be described as

j(t,pres;A,ω)=LL(v(x,r(t);pres)ṙ(t))dx,
(36)

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 vpress(x) = vp(x, r(t)) and the wall motion-driven flow (the solution of Equations (11)-(14)) as vpump(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

vpress(x,y=r(t))=κ(PL+PL2pres2coshλ(t)x/r(t)cosh(Lλ(t)/r(t))+PLPL2sinhλ(t)x/r(t)sinhLλ(t)/r(t)),
(37)

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

A(t)κμr(t)=121cos2(λ(t))sin(λ(t))λ(t)cos(λ(t)).
(38)

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 A0, and the only difference in Equation (37) and the asymptotic approximations is that one sets λ(t)3A(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

t=t̃/γ,pres=αr̄2p̃res,r(t)=r̃(t)r̄,
(39)

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

j̃press(t,pres)=r̃(t̃)λ(t̃)(ΞΥp̃res(t̃))tanh(λ(t̃)(t̃)),
(40)

where

r̃(t̃)=1+Θsin(Ωt̃),
(41)

and the nondimensional constants are defined as

Ξ=κ(PL+PL)γr̄,Υ=2αr̄κγ,Ā=κμr̄,
̄=L/r̄,Θ=Ar̄,Ω=ωγ,
(42)

and the time varying nondimensional parameters are defined as

A(t̃)=Ā/r̃(t̃),(t̃)=̄/r̃(t̃).
(43)

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 pres presented in Equation (1), we will begin by examining an intermediate result in which pres is a fixed value, pres(t) = Pres.

In the case when pres(t) = Pres, 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,

jj0=limTTT+2π/Ωj(t,pres)dt2πj0/Ω
(44)

where j0 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

j̃pressj̃press0=limTΩ2πTT+2π/Ωr̃(t̃)3/2tanh(λ(t̃)(t̃))tanh(λ̄̄)dt̃,
(45)

where j̃press0 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 Ā[1014,106]. These two constraints determine a maximal value for ̄ of 265, to ensure we remain within the radius of convergence of the formal power series for vpump. With these constraints in mind, we proceed by setting Ā=106, ̄=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.

FIG. 3.

We plot the percent change away from the non-pumping case of instantaneous flux as a function of the changing wall for a variety of values for Θ (a). In this case τ ∈ [2πn/Ω, 2π(n + 1)/Ω] where n ∈ ℤ. We average the change in flux over τ and plot the result as a function of Θ (b). Other parameters are fixed at Ā=106 and ̄=200.

FIG. 3.

We plot the percent change away from the non-pumping case of instantaneous flux as a function of the changing wall for a variety of values for Θ (a). In this case τ ∈ [2πn/Ω, 2π(n + 1)/Ω] where n ∈ ℤ. We average the change in flux over τ and plot the result as a function of Θ (b). Other parameters are fixed at Ā=106 and ̄=200.

Close modal

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

p̃̇res=p̃res+r̃(t̃)λ(t̃)(ΞΥp̃res(t̃))tanh(λ(t̃)(t̃))+ΘΩcos(Ωt̃)(t̃)(t̃)(vpump(x,1;A(t̃),(t̃))1)dx,
(46)

where vpump is the (nondimensional) solution for Equations (19)-(22) parametrized by A(t̃) and (t̃). 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

p̃res(t̃)=p̃res(0)exp0t̃g(s)+0t̃expst̃g(r)drf(s)ds
(47)

where

g(t̃)1+Υr̃(t̃)λ(t̃)tanh(λ(t̃)(t̃)),
(48)
f(t̃)Ξr̃(t̃)λ(t̃)tanh(λ(t̃)(t̃))+ΘΩcos(Ωt̃)(t̃)(t̃)(vpump(x,1;A(t̃),(t̃))1)dx.
(49)

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̃). We note that average value of f(t̃) 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̃) is small compared to g(t̃)1, i.e., ΘΩ/Υ = ϵ ≪ 1 (since Ξ = 0). We have determined in Sec. II that f(t̃) is absolutely bounded within its radius of convergence, and further note that g(t̃)>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̃), 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̃. 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, Ā=106, ̄=200, and vary |Θ| from 0.1 to 0.5, and Υ from 0 to 0.05. We then determine the long time average of p̃res(t̃),

p̃res(t̃)=Ω2πlimTTπ/ΩT+π/Ωp̃res(t̃)dt̃.
(50)

The reason we seek to determine the long time average of p̃res(t̃) is that it is the average metabolic material flux. Additionally it is proportional to the pressure driven flux. Thus if p̃res(t̃) 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 vpump 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 Θ.

FIG. 10.

Long time average reservoir pressure as a function of Υ for several values of Θ (a), Ω (b), ̄ (c), and Ā (d). By default Θ = 0.5, Ω = 10, Ā=106, ̄=200, and Ξ = 0.

FIG. 10.

Long time average reservoir pressure as a function of Υ for several values of Θ (a), Ω (b), ̄ (c), and Ā (d). By default Θ = 0.5, Ω = 10, Ā=106, ̄=200, and Ξ = 0.

Close modal
FIG. 4.

Long time average reservoir pressure as a function of Υ for several values of Θ. Ω = 10, Ā=106, ̄=200, and Ξ = 0.

FIG. 4.

Long time average reservoir pressure as a function of Υ for several values of Θ. Ω = 10, Ā=106, ̄=200, and Ξ = 0.

Close modal

We continue a sensitivity analysis in Appendix C 1 and find that Υ is also insensitive to changes in ̄ and Ā, 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 ̄, Ā, Ω, 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̃res has a stable steady state value given by

p̃ressslimtp̃res=Ξr̄λ(0)×tanh(λ(0)̄)1+Υtanh(λ(0)̄).
(51)

Similar to the above study in Section III B 1, we set Ā=106, ̄=200, Ω = 10, Θ = 0.5 and vary Υ; for each value of Υ, we set Ξ so that p̃resss=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 ̄ and Ā fixed, the average pressure is expected to converge to p̃resss=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̃resss1.

FIG. 5.

We fix p̃resss=1 and vary Υ setting Θ = 0.5, Ā=106, ̄=200, and Ω = 10.

FIG. 5.

We fix p̃resss=1 and vary Υ setting Θ = 0.5, Ā=106, ̄=200, and Ω = 10.

Close modal
TABLE I.

We examine the sensitivity of p̃resss with respect to Ω and Θ with p̃resss=1, Υ = 10−5, ̄=200, and Ā=106. When not being varied, Ω = 10 and Θ = 0.5.

Ω 10 0.1 0.01 
p̃res 0.9631 0.9623 0.9614 0.9614 
Θ 0.4 0.3 0.2 0.1 
p̃res 0.9791 0.9894 0.9956 0.9989 
Ω 10 0.1 0.01 
p̃res 0.9631 0.9623 0.9614 0.9614 
Θ 0.4 0.3 0.2 0.1 
p̃res 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

p̃res(t̃)=p̃res(0)exp0t̃g(s)+Ξ0t̃expst̃g(r)drh(Ωs)ds+O(Ω),
(52)

where

h(Ωt̃)Ξg(t̃)1Υ.
(53)

To study this solution, we first ignore the transient term by setting p̃res(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̃=2πτ/Ω 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̃. We find that as Ω → 0,

p̃res(2πτ/Ω)=02πτ/Ωesh(2πτΩs)ds02πτ/Ωesh(2πτ)ds=h(2πτ)=Ξg(2πτ/Ω)1Υ.
(54)

Figure 6 shows p̃res as a function of τ for several values of Θ, Ā, and ̄. 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.

FIG. 6.

We fix p̃resss=1 and vary τ over a full period of expansion and contraction. When not being varied, the parameters are fixed at Θ = 0.5, Ā=106, and ̄=200.

FIG. 6.

We fix p̃resss=1 and vary τ over a full period of expansion and contraction. When not being varied, the parameters are fixed at Θ = 0.5, Ā=106, and ̄=200.

Close modal

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 respiration4) or in a wave like pattern (renal filtration3). 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.

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.

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

Proposition.
If Equations (26)-(28) are closed with a normal wall velocity prescribed by annth-degree polynomial ofx2,
vn(x,±1)=±j=0najx2j,
then the pressurepnat the boundary will be an (n + 1)th degree polynomial ofx2,
pn(x,±1)=±j=0n+1bjx2j.

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 vn with j < n − 1 leads to a residual pressure that is a polynomial of degree at most (n − 1) of x2 at the boundary. Therefore we must demonstrate that

vn(x,±1)=±anx2(n1)
(A1)

leads to a residual pressure that is a polynomial of degree at most n of x2 at the boundary. We set an = 1 without loss of generality, due to the linearity of the equations.

Let ψn be the flow potential such that un = ∂yψn and vn = − ∂xψn. Due to the symmetry of the boundary conditions on the pressure and velocity at the boundaries, we assume that the flow profile un is odd in x and even in y, that vn is even in x and odd in y, and that pn is even in both x and y. Next, let us suppose there exists a finite power series such that

ψn=j=0(n+1)i=0(n+1j)cijx2i+1y2j+1.
(A2)

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

p(x,y)=Δyϕndx+f(y),
(A3)

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

ΔΔψn=0,
(A4)
xψn(x,±1)=±x2(n1),
(A5)
yψn(x,±1)=0.
(A6)

The only unknowns left in the system are the cij’s of which there are (n + 3)(n + 2)/2. Plugging ψn into the biharmonic equation and equating powers x and y results in (n + 1) n/2 equations described by

P4(2i+3)c(i+1)(j1)+2P2(2i+1)P22j+1cij+P4(2j+3)c(i1)(j+1)=0,
(A7)

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

j=0n+1i(2i+1)cij=δin,
(A8)
j=0n+1i(2j+1)cij=0
(A9)

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,

(2n+3)c(n+1)0=0,and
(A10)
c(n+1)0=0.
(A11)

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

13(2n+1)(2n+1)cn0cn1=01,
(A12)

which has solution

cn0=32(2n+1),cn1=12(2n+1).
(A13)

Proceeding recursively, we assume that we have found cij 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 cij for j ≥ 2, thereby reducing the (n + 2 − i) to two parameters (namely, ci0 and ci1). From here, we note that Equations (A8) and (A9) are reduced to

13(2i+1)(2i+1)ci0ci1+C1C2=00,
(A14)

where

C1=j=2n+1i(2j+1)cij,C2=j=2n+1i(2i+1)cij,
(A15)

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 cij (for 1 < j < n + 1 − i) by setting it as a function of c(i+1)(j−1) and c(i+2)(j−2),

cij=P4(2i+3)c(i+2)(j2)+2P2(2i+3)P22j1c(i+1)(j1)P4(2j+1).
(A16)

This procedure may then be repeated up to the case in which i = 1 and we have thus found a recursive algorithm for the solution to Equations (26)-(28) with normal velocity boundary condition given by Equation (A1). This completes the proof of Proposition 1.■

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

v(x,±r(t))=ṙ(t)/c,
(B1)

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

v(x,±1)=1/c.
(B2)

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

1c=v̄12v(x,1)dx.
(B3)

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

σsi1|v̄|12(v(x,1)v̄)2dx.
(B4)

σ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 A2 is fixed. To see this concretely, we overlay plots of σsi as a function of A2 for values of A[106,3.16×103] in Figure 8.

FIG. 7.

The variation σsi away from suction/injection boundary conditions. Note that for all values of the nondimensional permeability, A, the assumption suction/injection boundary conditions breaks down as the ratio between channel length and width, ℓ, is sufficiently large. Additionally, we display the level set curve on which the variation is 10% (white dot-dashed line) and 20% (white dashed line).

FIG. 7.

The variation σsi away from suction/injection boundary conditions. Note that for all values of the nondimensional permeability, A, the assumption suction/injection boundary conditions breaks down as the ratio between channel length and width, ℓ, is sufficiently large. Additionally, we display the level set curve on which the variation is 10% (white dot-dashed line) and 20% (white dashed line).

Close modal
FIG. 8.

Variation σsi as a function of A2 is in roughly invariant for A[106,3.16×103]. The profile is nearly invariant for these small values of A.

FIG. 8.

Variation σsi as a function of A2 is in roughly invariant for A[106,3.16×103]. The profile is nearly invariant for these small values of A.

Close modal

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

x̃=x/r(t),ỹ=y/r(t),
(B5)
ũ=u/ṙ(t),ṽ=v/ṙ(t),
(B6)
p̃=pκ/Aṙ(t).
(B7)

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

p(,±1)=0,p(,±1)=0,
(B8)
p=Δu,u=0,
(B9)
u(x,±1)=0,
(B10)
v(x,±1)=±(Ap(x,±1)+1).
(B11)

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

u(x,y)320x1+y2×10+A1152+5x25y2,
(B12)
v(x,y)120y[103+y2+3A{523+y25x23+y2+1+y22}].
(B13)

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

σp(w)supx1,x2[0,]infβR11(w(x1,y)βw(x2,y))2dy11w(x1,y)2dy.
(B14)

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 β=maxy(|w(x1,y)/w(x2,y)|) and consider the limit as x1 → ℓ and x2 → 0. Our choice for β ensures that computed value will be an upper bound for σp; our choice for x1 and x2 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.

FIG. 9.

We plot max(σp(u), σp(v)) for the asymptotic expansion presented in Equations (B12) and (B13) (a). We also plot max(σp(u), σp(v)) for the solutions that have converged numerically with an error tolerance of 10−4 (b).

FIG. 9.

We plot max(σp(u), σp(v)) for the asymptotic expansion presented in Equations (B12) and (B13) (a). We also plot max(σp(u), σp(v)) for the solutions that have converged numerically with an error tolerance of 10−4 (b).

Close modal

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 p̃res as a function of Υ for Ω ∈ [10, 100], ̄[100,200], and Ā[108,106] (see Figures 10(b)-10(d), respectively); additionally, Θ = 0.5, Ω = 10, Ā=106, and ̄=200. Our results indicate that the optimal Υ is insensitive with respect to Θ and Ā, changes with respect to ̄, 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, Ā=107, ̄=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.

FIG. 11.

We plot the long time average reservoir pressure as a function of Υ for several values of Θ (a), Ω (b), ̄ (c), and Ā (d). When the variables are not being varied, we fix Θ = 0.5, Ω = 0.1, Ā=107, ̄=500, and Ξ = 0.

FIG. 11.

We plot the long time average reservoir pressure as a function of Υ for several values of Θ (a), Ω (b), ̄ (c), and Ā (d). When the variables are not being varied, we fix Θ = 0.5, Ω = 0.1, Ā=107, ̄=500, and Ξ = 0.

Close modal

a. The dependency of p̃res on Ā, ̄, 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 Ā, ̄, and Ω is shown in Figure 12. By default Ā=106, ̄=200, Ω = 10, Υ = 0.02, Θ = 0.3, and Ξ = 0. Our results indicate that the mean reservoir pressure increases with Ω, Ā, and ̄. 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).

FIG. 12.

We plot the long time average reservoir pressure as a function of Ā, ̄, and Ω. When not being viewed as an independent variable, Ā=106, ̄=200, and Ω = 10. We also fix Υ = 0.02, Θ = 0.3, and Ξ = 0.

FIG. 12.

We plot the long time average reservoir pressure as a function of Ā, ̄, and Ω. When not being viewed as an independent variable, Ā=106, ̄=200, and Ω = 10. We also fix Υ = 0.02, Θ = 0.3, and Ξ = 0.

Close modal

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, Ā=106, ̄=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 p̃res0 when Υ → 0 or Υ → ∞. Thus, without changing scales, we propose three hypotheses as to why we observe enhanced material extraction:

  1. Pumping generates material extraction because of the skewness of the motion driven flux (f(t̃)).

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

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

FIG. 13.

We plot the coefficients on the differential equation for the dynamic reservoir model (Equation (1)) with Ā=106, ̄=200, and Ω = 10, Υ = 0.01, Θ = 0.5, and Ξ = 0. f and g are defined in Equations (49) and (48), respectively.

FIG. 13.

We plot the coefficients on the differential equation for the dynamic reservoir model (Equation (1)) with Ā=106, ̄=200, and Ω = 10, Υ = 0.01, Θ = 0.5, and Ξ = 0. f and g are defined in Equations (49) and (48), respectively.

Close modal

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

y=g(t)y+F(t;s),
(C1)

where

F(t;s)=c1H(t;s)H(t;s),
(C2)

where c1 is some constant,

H(t;s)=tanhehΩt2πΩt+π2π2tanh[eh/2]tanhehΩt2πΩt+π2π2tanh[eh/2],
(C3)

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 c1 = 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.

FIG. 14.

We plot the (unitless) amount of material extraction when considering Equation (C1). Ā=106, ̄=200, and Ω = 10, Υ = 0.01, Θ = 0.5, Ξ = 0, and c1 = 50.

FIG. 14.

We plot the (unitless) amount of material extraction when considering Equation (C1). Ā=106, ̄=200, and Ω = 10, Υ = 0.01, Θ = 0.5, Ξ = 0, and c1 = 50.

Close modal

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

y=G(t;r)y+F(t;0),
(C4)

where

G(t;r)=c2+h(t;r)tt+2π/Ωh(τ;r)dτ,
(C5)

c2 is some constant,

h(t;r)=log(1+rsin(Ωt))log(1+r)log(1r),
(C6)

and r ∈ (0, 1) is a parameter that controls the symmetry of G about its mean c2. 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 c2 = 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).

FIG. 15.

We plot the (unitless) amount of material extraction when considering Equation (C4). Ā=106, ̄=200, and Ω = 10, Υ = 0.01, Θ = 0.5, Ξ = 0, c1 = 50, and c2 = 1.

FIG. 15.

We plot the (unitless) amount of material extraction when considering Equation (C4). Ā=106, ̄=200, and Ω = 10, Υ = 0.01, Θ = 0.5, Ξ = 0, c1 = 50, and c2 = 1.

Close modal

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

y=(2+Υsin(Ωt))y+sin(Ωt+ϕ).
(C7)

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.

FIG. 16.

We plot the (unitless and normalized) amount of material extraction when considering Equation (C7); Ω = {0.1, 1, 10}, and Υ = 1 (a). We test the idea of the phase shift on the physical model for Ā=106, ̄=200, Ω = 10, Υ = 0.01, Θ = 0.5, Ξ = 0 (b).

FIG. 16.

We plot the (unitless and normalized) amount of material extraction when considering Equation (C7); Ω = {0.1, 1, 10}, and Υ = 1 (a). We test the idea of the phase shift on the physical model for Ā=106, ̄=200, Ω = 10, Υ = 0.01, Θ = 0.5, Ξ = 0 (b).

Close modal

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

y(t)=g(t+ψ)+f(t),
(C8)

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.

1.
T.
Weis-Fogh
, “
Respiration and tracheal ventilation in locusts and other flying insects
,”
J. Exp. Biol.
47
,
561
587
(
1967
).
2.
Y.
Komai
, “
Augmented respiration in a flying insect
,”
J. Exp. Biol.
201
,
2359
2366
(
1998
).
3.
B.
Schmidt-Nielsen
and
B.
Schmidt-Nielsen
, “
On the function of the mammalian renal papilla and the peristalsis of the surrounding pelvis
,”
Acta Physiol.
202
,
379
385
(
2011
).
4.
Y.
Aboelkassem
,
A. E.
Staples
, and
J. J.
Socha
, “
Microscale flow pumping inspired by rhythmic tracheal compressions in insects
,”
Proc. ASME
471
479
(
2011
).
5.
Y.
Aboelkassem
and
A. E.
Staples
, “
Flow transport in a microchannel induced by moving wall contractions: A novel micropumping mechanism
,”
Acta Mech.
223
,
463
480
(
2012
).
6.
G.
Hosseini
,
J.
Williams
,
E.
Avital
,
A.
Munjiza
,
X.
Don
, and
J.
Green
, “
Computational simulation of the urinary system
,” in
Proceedings of the world congress on Engineering and Computer Science
(
International Association of Engineers
,
2012
), Vol.
2201
, Issue 1, pp. 710–714.
7.
S. A.
Regirer
, “
On the approximate theory of the flow of a viscous incompressible liquid in a tube with permeable walls
,”
Zh. Tekh. Fiz.
30
,
639
643
(
1960
).
8.
L. S.
Galowin
,
L. S.
Fletcher
, and
M. J.
DeSantis
, “
Investigation of laminar flow in a porous pipe with variable wall suction
,”
AIAA J.
12
,
1585
1589
(
1974
).
9.
J.
Granger
,
J.
Dodds
, and
N.
Midoux
, “
Laminar flow in channels with porous walls
,”
Chem. Eng. J.
42
,
193
204
(
1989
).
10.
S. K.
Karode
, “
Laminar flow in channels with porous walls, revisited
,”
J. Membr. Sci.
191
,
237
241
(
2001
).
11.
P.
Haldenwang
, “
Laminar flow in a two-dimensional plane channel with local pressure-dependent crossflow
,”
J. Flu. Mech.
593
,
463
473
(
2007
).
12.
P.
Haldenwang
and
P.
Guichardon
, “
Pressure runaway in a 2d plane channel with permeable walls submitted to pressure-dependent suction
,”
Eur. J. Mech. B/Fluids
30
,
177
183
(
2011
).
13.
N.
Tilton
,
D.
Martinand
,
E.
Serre
, and
R.
Lueptow
, “
Incorporating Darcys law for pure solvent flow through porous tubes: Asymptotic solution and numerical simulations
,”
AIChE J.
58
,
230
244
(
2012
).
14.
B.
Bernales
and
P.
Haldenwang
, “
Laminar flow analysis in a pipe with locally pressure-dependent leakage through the wall
,”
Eur. J. Mech. B/Fluids
43
,
100
109
(
2014
).
15.
G.
Herschlag
,
J.-G.
Liu
, and
A. T.
Layton
, “
An exact solution for stokes flow in a channel with arbitrarily large wall permeability
,”
SIAM J. Appl. Math.
75
,
2246
2267
(
2015
).
16.
A. S.
Berman
, “
Laminar flow in channels with porous walls
,”
J. Appl. Phys.
24
,
1232
1235
(
1953
).
17.
S. W.
Yuan
and
A. B.
Finkelstein
, “
Laminar pipe flow with injection and suction through a porous wall
,”
Trans ASME
78
,
719
724
(
1956
).
18.
R. M.
Terrill
, “
Laminar flow in a uniformly porous channel
,”
Aeronaut. Q.
15
,
1033
1038
(
1964
).
19.
R.
Terrill
, “
An exact solution for flow in a porous pipe
,”
J. Appl. Math. Phys.
33
,
547
552
(
1982
).
20.
J.
Majdalani
,
C.
Zhou
, and
C.
Dawson
, “
Two-dimensional viscous flow between slowly expanding or contracting walls with weak permeability
,”
J. Biomech.
35
,
1399
1403
(
2002
).
21.
E.
Dauenhauer
and
J.
Majdalani
, “
Exact self-similarity solution of the Navier-Stokes equations for a porous channel with orthogonally moving walls
,”
Phys. Fluids
15
,
1485
(
2003
).
22.
S.
Asghar
,
M.
Mushtaq
, and
A.
Kara
, “
Exact solutions using symmetry methods and conservation laws for the viscous flow through expanding-contracting channels
,”
Appl. Math. Modell.
32
,
2936
2940
(
2008
).
23.
S.
Asghar
,
M.
Mushtaq
, and
T.
Hayat
, “
Flow in a slowly deforming channel with weak permeability: An analytical approach
,”
Nonlinear Anal.: Real World Appl.
11
,
555
561
(
2010
).
24.
S.
Mohyud-Din
,
A.
Yildirim
, and
S. A.
Sezer
, “
Analytical approach to a slowly deforming channel flow with weak permeability
,”
Z. Naturforsch. A
65
,
299
310
(
2010
).
25.
S.
Xin-Hui
,
Z.
Lian-Cun
,
Z.
Xin-Xin
,
S.
Xin-Yi
, and
Y.
Jian-Hong
, “
Flow of a viscoelastic fluid through a porous channel with expanding or contracting walls
,”
Chin. Phys. Lett.
28
,
044702
(
2011
).
26.
M. M.
Rashidia
and
S. M.
Sadri
, “
New analytical solution of two-dimensional viscous flow in a rectangular domain bounded by two moving porous walls
,”
Int. J. Comput. Methods Eng. Sci. Mech.
12
,
26
33
(
2011
).
27.
M.
Azimi
,
A.
Hedesh
, and
S.
Karimian
, “
Flow modeling in a porous cylinder with regressing walls using semi analytical approach
,”
Mech. Mech. Eng.
18
,
77
84
(
2014
).
28.
Sushila
,
J.
Singh
, and
Y.
Shishodia
, “
A reliable approach for two-dimensional viscous flow between slowly expanding or contracting walls with weak permeability using sumudu transform
,”
Ain Shams Eng. J.
5
,
237
242
(
2014
).
29.
Wolfram Research, Inc., Mathematica, version 8.0, Champaign, IL, 2010.