Transient flows in active porous media

TL;DR

This paper models the coupled advective and diffusive transport in a channel with stimulus-responsive permeability, revealing damped oscillations in steady state and discussing experimental observability.

Contribution

It provides an exact analytical solution for coupled transport in a stimulus-responsive porous medium, highlighting oscillatory behavior in steady state.

Findings

Existence of damped oscillations in flow reaching equilibrium.

Exact solution to the coupled transport problem.

Feasibility of experimental observation discussed.

Abstract

Stimuli-responsive materials that modify their shape in response to changes in environmental conditions -- such as solute concentration, temperature, pH, and stress -- are widespread in nature and technology. Applications include micro- and nanoporous materials used in filtration and flow control. The physiochemical mechanisms that induce internal volume modifications have been widely studies. The coupling between induced volume changes and solute transport through porous materials, however, is not well understood. Here, we consider advective and diffusive transport through a small channel linking two large reservoirs. A section of stimulus-responsive material regulates the channel permeability, which is a function of the local solute concentration. We derive an exact solution to the coupled transport problem and demonstrate the existence of a flow regime in which the steady state is reached via a damped oscillation around the equilibrium concentration value. Finally, the feasibility of an experimental observation of the phenomena is discussed. Please note that this version of the paper has not been formally peer reviewed, revised or accepted by a journal.

Equations

$${\partial }_{T}C+\mathbf{v}\cdot \nabla C=D{\nabla }^{2}C,$$

$${\partial }_{T}C+\frac{{\kappa }_0\Delta p}{\eta L}\frac{C(C_0)}{C_0}{\partial }_{X}C=D{\partial }_X^{2}C.$$

$$C(0)=0,C(L)=c_0.$$

$$x=X/L,c=C/C_0,\text{and}t=T(D/L^2)$$

$${\partial }_{t}c+\text{Pe}c(x_0){\partial }_{x}c={\partial }_x^{2}c,$$

$$c(0)=0,c(1)=1.$$

$$c(x,0)=0.$$

$$\text{Pe}\gamma {\partial }_x{c}_s(x)={\partial }_x^2{c}_s(x),$$

$$c_s(x)=\frac{e^{\text{Pe}\gamma x}-1}{e^{\text{Pe}\gamma }-1}.$$

$$\gamma =\frac{e^{\text{Pe}\gamma x_0}-1}{e^{\text{Pe}\gamma }-1}.$$

$$c(x,t)=x+\frac{2}{\pi }\sum _n\frac{(-1{)}^n}{n}{\exp }^{-n^2{\pi }^{2}t}\sin (n\pi x).$$

$$c(x_{\text{obs}},t)-c_s(x_{\text{obs}})\propto e^{-(k_r+i{k}_i)t}$$

$$c(x,t)=c_s(x)-c_1(x,t),$$

$$c_1(x,t)=e^{-kt}g(x),$$

$$g^{\prime \prime }(x)-\text{Pe}\gamma g^{\prime }(x)+kg(x)=\frac{{\text{Pe}}^2\gamma }{e^{\text{Pe}\gamma }-1}{e}^{\text{Pe}\gamma x}g(x_0),$$

$$g(0)=0,g(1)=0,g(x_0)=1$$

$$\textrm{Pe}^{*}&=\textrm{Pe}\,\gamma,\\ \lambda=\frac{1}{k}\frac{\textrm{Pe}^{*2}}{\gamma(e^{\textrm{Pe}^{*}}-1)},\\ \zeta=\frac{1}{2}\textrm{Pe}^{*}\sqrt{\frac{4}{\lambda\gamma(e^{\textrm{Pe}^{*}}-1)}-1}.$$

$$g(x)=\lambda e^{{\text{Pe}}^{\ast }x}+e^{\frac{{\text{Pe}}^{\ast }x}{2}}(A\cos \zeta x+B\sin \zeta x)$$

$$g(0)=0,g(1)=0,g(x_0)=1,$$

$$\lambda+A&=0\\ \lambda e^{\textrm{Pe}^{*}}+e^{\frac{\textrm{Pe}^{*}}{2}}(A\cos\zeta+B\sin\zeta)=0\\ \lambda e^{\textrm{Pe}^{*}\,x_{0}}+e^{\frac{\textrm{Pe}^{*}\,x_{0}}{2}}(A\cos\zeta x_{0}+B\sin\zeta x_{0})=1\qquad$$

$$\lambda e^{\frac{{\text{Pe}}^{\ast }{x}_0}{2}}\left[e^{\frac{{\text{Pe}}^{\ast }{x}_0}{2}}-\cos \zeta x_0+\frac{\sin \zeta x_0}{\sin \zeta }\left(\cos \zeta -e^{\frac{{\text{Pe}}^{\ast }}{2}}\right)\right]=1.$$

$$g(x)=\lambda\bigg{[}e^{\textrm{Pe}^{*}\,x}+e^{\frac{\textrm{Pe}^{*}\,x}{2}}\bigg{(}-\cos\zeta x+\\ +\bigg{(}\cot\zeta-\frac{e^{\frac{\textrm{Pe}^{*}}{2}}}{\sin\zeta}\bigg{)}\sin\zeta x\bigg{)}\bigg{]},$$

$${\text{Pe}}^{\ast }=\frac{v}{v_0}\text{Pe}=\gamma \text{Pe}$$

$$k(Pe,x_0{)}^{1/2}=\pi +2\text{Pe}\frac{\sin (\pi x_0)}{{\pi }^2}+\mathcal{O}({\text{Pe}}^2).$$

Full text

I Introduction

Fluid flow and convective solute transport in porous media and confined channel geometries are ubiquitous in nature and technology. Interesting phenomena arise when channels walls and solid structures are themselves active; for instance, when the presence of solutes influences the channel geometry and hence permeability to fluid flow. Man-made examples include sensing and actuation in microfluidic systems using stimuli-responsive hydrogels Koetting et al. (2015). Responsive biomaterials are found, for example, in the phloem and xylem vascular systems of plants, where neighboring cells are separated by planar membranes covered with pores which respond to changes in concentration of chemical signals Zwieniecki et al. (2001); Mullendore et al. (2010). The stimuli that induce changes in these synthetic and natural materials have been widely studies. However, the coupling between induced volume changes and advective solute transport in porous materials, however, are not well understood.

In this paper we investigate the transient nature of advective transport in active porous media. We study a one-dimensional system where the advective solute transport speed is coupled to the concentration field. Numerical investigation of the model reveals the existence of a flow regime in which the steady state is reached via a damped oscillation around the equilibrium concentration value. We derive an exact solution using perturbation theory and show that the flow dynamics depends primarily on the ratio of advective to diffusive transport timescales (the Peclet number; Pe). Above a critical Pe-value, damped oscillations occur in both the velocity and concentration fields. Finally, we propose an experimental design to test the theoretical predictions.

II Flows in active porous media

Stimulus-responsive hydrogels have been a topic of extensive research the past decades Koetting et al. (2015). Their ability to modify their internal structure based on external stimuli allows for dynamic control over flows in biological Zwieniecki et al. (2001); Mullendore et al. (2010) or man-made systems Koetting et al. (2015); Huglin (1989) (Fig. 1). Responsive hydrogels, i.e. hydrophilic polymers embedded and crosslinked into hydrophilic structures Huglin (1989), can respond to a broad range of stimuli, e.g. pH Gupta et al. (2002); Sharpe et al. (2014); Pep (2000); Park et al. (1998); Chan et al. (2008); Chan and Neufeld (2009); Kim et al. (2000); Qu et al. (2000); Schmaljohann (2006), temperature Klouda and Mikos (2008); Purushotham and Ramanujan (2010), individual molecules (chemically driven) Wu et al. (2011); Bernfeld and Wan (1963); Horbett et al. (1984); Klumb and Horbett (1992), shear stress Hoffman (2013); Qiu and Park (2001); Wang et al. (2008); Bell et al. (2006); Guvendiren et al. (2012); Gutowska et al. (2001) etc., that trigger a change of material properties. In pH induced responses, hydrogel swelling/deswelling occurs when polymers are ionized by the dynamically changing environmental pH Gupta et al. (2002). Hence, the charge buildup results in an electrostatic force generation within the hydrogel that ultimately leads to absorbance or expulsion of water Sharpe et al. (2014); Pep (2000). Other workers have investigated temperature dependent hydrogels utilizing the critical solubility temperature with applications in drug delivery Klouda and Mikos (2008) and tissue engineering Purushotham and Ramanujan (2010). Hydrogels also exhibit responsive behaviour to chemical stimuli such as glucose by entrapping glucose oxidase enzymes in the hydrogel structure Wu et al. (2011). Another group of stimulus responsive hydrogels are known to respond to mechanical stress. Two subgroups that emerge are materials with shear thinning or shear thickening behaviour due to the viscoelastic nature of systems comprised of polymers, an intermediate material state at the interface between liquids and solids Hoffman (2013); Qiu and Park (2001). Applications of shear stress responsive hydrogels include, among others, drug delivery and wound repair Wang et al. (2008); Bell et al. (2006); Guvendiren et al. (2012); Gutowska et al. (2001).

In summary, the physiochemical factors that induce volume chances in stimulus-responsive materials are well understood. By contrast, less is known about the coupling between fluid flow, solute advection and stimulus response in these systems.

III Model

To elucidate the transient behavior of flow in active porous media, we consider flow in a small channel of constant cross section aligned with the horizontal $X$-axis linking two large reservoirs (Fig. 2). The channel has length $L$, width $w$ and height $h$, and a short section of active porous media is located at $X=X_{0}<L$. The right reservoir ($X=L$) is kept at constant concentration $X_{0}$, while the left reservoir ($X=0$) contains no solute. This drives a diffusive flux $-D\nabla C$ in the channel where $D$ is the diffusion coefficient. The right reservoir ($X=L$) is kept at constant pressure $p_{0}$, while the left reservoir ($x=0$) is at a higher pressure $p_{0}+\Delta p$. We assume the advective flow speed in the channel $v$ follows Darcy’s law, $v\sim\kappa\Delta p/(\eta L)$, where $\kappa$ is the channel permeability and $\eta$ is the viscosity. To model the active porous media, we assume that the dependence of the channel permeability $k$ on solute concentration $c$ can be expressed as $\kappa=\kappa_{0}C(X_{0})/C_{0}$. The hydraulic conductivity is proportional to the concentration at location the active porous media, $X_{0}$, such that high solute concentration evokes deswelling of the post valves while low concentration to an increase of post valve volume (Fig. 2).

The transport of solutes in the channel is governed by the advection-diffusion equation

[TABLE]

where $T$ is time, $\mathbf{v}$ is the velocity field and $D$ is the diffusion coefficient. With the aforementioned assumptions, this reduces to a one-dimensional equation for the concentration $C(X,T)$ in the channel

[TABLE]

The boundary conditions are

[TABLE]

For convenience we introduce the non-dimensional variables

[TABLE]

The dimensionless governing equation is

[TABLE]

where $x_{0}=X_{0}/L$ and we have introduced the dimensionless Peclet number $\textrm{Pe}=v_{0}L/D$. Here, $v_{0}=\kappa_{0}\Delta p/(\eta L)$ is the maximum reference velocity. The Peclet number characterizes the relative contribution from advective and diffusive transport. The boundary conditions in Eq. (3) become

[TABLE]

In the following, we consider the initial condition corresponding to an empty channel:

[TABLE]

and study the transient dynamics of the system. Before proceeding, however, we briefly discuss the steady-state solution $c_{s}(x)$ to Eq. (5) and the system behavior when $\textrm{Pe}\ll 1$.

III.1 Steady-state solution

When $\partial_{t}c=0$, Eq. (5) reduces to

[TABLE]

where we have introduced the parameter $\gamma=c_{s}(x_{0})$, the steady-state concentration at $x_{0}$. The solution to Eq. (8) with boundary conditions (6) is

[TABLE]

The parameter $\gamma=c_{s}(x_{0})$ can be determined as a function of the system parameters ($Pe$ and $x_{0}$) by solving the trancendental equation

[TABLE]

When $\textrm{Pe}\ll 1$, we find that $c(x)=x$ and $\gamma=x_{0}$. Taking the limit $\textrm{Pe}\gg 1$ leads to $\log\gamma/\gamma=\textrm{Pe}(x_{0}-1)$

III.2 Solution in a diffusion dominated system

When diffusion dominates, the Peclet number $\textrm{Pe}\ll 1$. In that limit, the solution to Eqns. (5)-(6) is

[TABLE]

Equilibrium is approached exponentially on the timescale set by the slowest mode ($n=1$). The system is within $5\%$ of the steady state solution when $t\simeq 3\pi^{-2}\simeq 0.3$.

IV Results

IV.1 Numerical Simulation

In order to reveal the transient nature of flow in active porous media (Fig. 2), we ran simulations of Eqns. (5)-(7) for a range of values for Pe and $x_{0}$. For relatively low Peclet numbers – corresponding to a diffusion-dominated system – the steady state is reached asymptotically with non-dimensional relaxation time $\tau\sim 0.3$. The behavior of the system is thus in accord with a purely diffusive process (Eq. (11)), where equilibrium is approached exponentially on a similar timescale.

By contrast, for values of the Peclet number Pe above unity, the characteristics of the system changes in two respects. First, the steady state is reached on a timescale which decreases with increasing Pe. Second, for large Pe the approach to equilibrium follows a damped oscillation (Fig. 3(c)), indicating a qualitative deviation from the asymptotic approach to equilibrium in Eq. (11) and Fig. 3(a).

To further elucidate the characteristics of the oscillations, we studied the temporal evolution of a small disturbance to the steady state. We thus added a weak gaussian perturbation to the steady-state solution (Eq. (9)) at an arbitrary position $x_{p}$ within the domain, and studied the approach to equilibrium. After the initial perturbation had decayed, we observed an approximately decaying harmonic time dependence of the disturbance at the position $x_{\text{obs}}$, i.e.

[TABLE]

where $c_{s}$ is the steady state concentration given in Eq. (9). In Eq. (12), $k_{r}$ and $k_{i}$ are the real and imaginary parts of the complex wavenumber $k$, corresponding to decay time $\sim k_{r}^{-1}$ and oscillation period $\sim 2\pi k_{i}^{-1}$. We thus extracted $k_{r}$ and $k_{i}$ from the numerical simulations by curve fitting using. Neither the position of the perturbation $x_{p}$ nor the observation location $x_{\text{obs}}$ appeared to influence the magnitude of the wavenumber $k$ significantly. However, we chose the parameters to avoid overlap between the position of the active porous material $x_{0}$ and $x_{p}$ and $x_{\text{obs}}$. Finally, we found that while oscillations are present when the position is to the right of the channel centerline ($x_{0}>1/2$), they decay rapidly and a sum of at least two decaying exponentials are necessary to provide a satisfactory curve fit. In the following, we thus restrict ourselves to the case $x_{0}<1/2$.

Having extracted the wavenumber $k=k_{r}+ik_{i}$ from the numerical simulations, we studied $k$’s dependence on the relative importance of advection and diffusion (Fig. 4). When the Peclet number is relatively small, we found $k_{r}\simeq 10$ and $k_{i}=0$, in accord with Eq. (11), which predicts $k_{r}=1^{2}\pi^{2}\simeq 9.87$ and $k_{i}=0$. The simulations further revealed that the onset of oscillations occurs at a critical value of the Peclet number $\textrm{Pe}_{c}$. For the case $x_{0}=1/4$ shown in Fig. 4, $\textrm{Pe}_{c}\simeq 20$. Note that the magnitude of the critical $\textrm{Pe}_{c}$ varies depending on the location $x_{0}$ of the active porous media in the channel (see also Fig. 6).

The physical mechanism that triggers the onset of oscillations can be interpreted as follows: when the system is perturbed away from the steady state $c_{s}(x)$, the concentration at $x_{0}$, i.e. $c(x_{0})$, will shift either up or down as solute is transported by a convention through the domain. This directly influences the advective flow speed, which is proportional to the local concentration at that point, $v\propto c(x_{0})/c_{0}$. Diffusion will counteract this process, eventually returning the system to the steady state $c_{s}(x)$. However, if the advective transport is sufficiently strong, advection can push the system into a state in where the concentration gradients become so great that the concentration $c(x_{0})$ overshoots it’s equilibrium value as diffusion counteracts advection. The process repeats itself – with a progressively smaller amplitude – until the steady state is restored.

IV.2 Analytic Solution

To rationalize the observed onset of oscillations at high Peclet-numbers (Fig. 3 and 4) and their dependence on the system parameters, we proceed to consider the evolution of the perturbed system. Considering a small deviation from the steady state $c_{s}(x)$, we write

[TABLE]

where we assume the perturbation $c_{1}\ll c_{s}$. We further assume that the perturbation has a harmonic time dependence

[TABLE]

where $k=k_{r}+ik_{i}$ is the complex wavenumber and $k_{r}>0$ and $g(x)$ is an unknown function of $x$. Substitution of Eq. (14)-(13) into Eq. (5) leads to a spatial equation for $g(x)$

[TABLE]

where prime denotes derivative with respect to $x$. The boundary conditions are

[TABLE]

where we have eliminated quadratic terms in $g$. Note that $g(x_{0})$ is an arbitrary constant that defines the strength of the perturbation, chosen here as unity.

Equation (15) is solved following the method of Pedley and Fischbarg (1978), who analyzed a similar problem related to transient flows near osmotic membranes. A particular solution to the inhomogeneous equation is $g_{i}(x)=\lambda\exp\left(\textrm{Pe}^{*}x\right)$, while the homogeneous solution is $g_{h}=\exp\left(\textrm{Pe}^{*}x/2\right)(A\cos\zeta x+B\sin\zeta x)$. Here, we have introduced the parameters

[TABLE]

The complete solution to (15) is

[TABLE]

To determine the constants $A$ and $B$ we apply the boundary conditions

[TABLE]

which after substitution become

[TABLE]

By eliminating $A$ and $B$ we find an eigenvalue equation for the wavenumber $k$:

[TABLE]

To test the validity of our solution, we compared the predictions from Eq. (21) with numerical data. For a given set of parameters $(\textrm{Pe},x_{0})$ we thus determined the solution to Eq. (21) with the smallest real part of $k$, corresponding to the slowest decaying mode. The solutions to Eq. (21) are in good agreement with the numerically extracted eigenvalues (Fig. 4). The eigenfunction is

[TABLE]

shown in Fig.  5. We note that the spatial eigenfunctions in Eq. (22) are consistent with Eq. (11) when Pe is relatively small.

IV.3 Critical Pe for onset of oscillations

To elucidate the conditions under which damped oscillations occur in our system, we extracted a phase diagram (Fig. 6) from the eigenvalue equation (21). Oscillations in the mode associated with the smallest real eigenvalue can occur for values for the Peclet number at or above $\textrm{Pe}\simeq 18$, depending on the position of the active porous media $x_{0}$. This suggests that advection should be nearly twenty times stronger than diffusion to obtain oscillations. However, because of the coupling between the permeability of the porous media $\kappa=\kappa_{0}c(x_{0})/c_{0}=\kappa\gamma$ and concentration $c(x_{0})$, we can write for the flow speed $v=\kappa\Delta p/(\eta L)=\gamma v_{0}$. This implies that the physically relevant Peclet number is $\textrm{Pe}^{*}$, given by

[TABLE]

Replotting the phase diagram using the rescaled Peclet number $\textrm{Pe}^{*}$ reveals that the onset of oscillations occur when advection is $2$ to $4$ times stronger than diffusion.

Oscillations are found in the numerical simulations for $x_{0}>1/2$. As noted earlier, however, they decay rapidly and a sum of at least two decaying exponentials are necessary to provide a satisfactory curve fit. By extracting the three smallest roots of Eq. (21), we found that for $x_{0}>1/2$ the oscillations are no longer associated with mode with the smallest real eigenvalue. Oscillations are found, however, in higher-order solutions to Eq. (21), an observation which provides a qualitative rational for the numerical results.

IV.4 Small-Pe expansion

We end this section by deriving an analytical expression for the solution to Eq. (21) for small Pe. Taking the limit $\textrm{Pe}\ll 1$ in Eq. (20) leads to $\lambda=\textrm{Pe}/k$ and $\zeta=\sqrt{k}$. Interting this into Eq. (21) and assuming that we can write the eigenvalue as a power-series in Pe: $k^{1/2}=a_{0}+\textrm{Pe}\,a_{1}+\mathcal{O}(\textrm{Pe}^{2})$ with $a_{0}=\pi$ gives an analytical expression for the eigenvalue $k$ at low-Pe

[TABLE]

The approximate expression in Eq. (24) is in reasonable accord with the solution to Eq. (21) for $Pe\lesssim 10$ (Fig. 4).

V Discussion and conclusion

A relatively complete picture of the factors that influence transient flows in active porous media has emerged. First and foremost, we have demonstrated the existence of damped oscillations in the flow within a channel linking two reservoirs (Fig. 2). The oscillations in solute concentration and liquid flow speed is the result of a coupling between solute concentration $c$ and the permeability of the channel, controlled by the swelling and shrinking of a stimulus-responsive material located at the position $x_{0}$ in the channel. Damped oscillations occur when the advective transport $\sim vc$ is sufficiently great to overcome diffusive transport $\sim Dc/L$, i.e. when the Peclet number $\textrm{Pe}^{*}=vL/D$ is greater than a critical value, which varies in the range from $\textrm{Pe}^{*}\simeq 2$ to $\textrm{Pe}^{*}\simeq 4$, depending on the location of the active porous material $x_{0}$ (Fig. 6).

To observe the damped oscillations in a laboratory setting, we an experiment based on the system in Fig. 2. In a channel of length $L=1$ mm, the characteristic diffusive time is $t_{d}=L^{2}/D=2\times 10^{3}$ s, where we have used the diffusion coefficient $D=5\times 10^{-10}$ m2/s of the dye carboxyflourescein Carroll et al. (2014). With a stimulus-responsive hydrogel valve located at $x_{0}=0.25$ mm operated at a moderate Peclet number of $\textrm{Pe}^{*}=4$ ($\textrm{Pe}=100$), we find $k_{r}\sim 30$ and $k_{i}\sim 20$. This corresponds to an oscillation period of $T_{\text{osc}}=2\pi k_{i}^{-1}t_{d}\simeq 630$ s while the decay time is $\tau=3t_{d}k_{r}^{-1}\simeq 200$ s. At least the first half-period should be observable for this choice of parameters. For a slightly lower forcing $\textrm{Pe}^{*}=2$ ($\textrm{Pe}\simeq 30$), we find $k_{r}\sim 30$ and $k_{i}\sim 10$ corresponding to oscillation period $T_{\text{osc}}\simeq 1300$ s and decay time $\tau\simeq 200$ s. For water flowing in a channel of width $w=300$ $\mathrm{\mu m}$ and height $h=100$ $\mathrm{\mu m}$, these flows would require pressure differentials of $\Delta p=0.077$ Pa and $\Delta p=0.023$ Pa respectively Bruus (2007). These pressure differentials, in the absence of post valves would generate flow velocities $v_{0}=\textrm{Pe}D/L=50$ $\mathrm{\mu m/s}$ and $v_{0}=15$ $\mathrm{\mu m/s}$ for $\textrm{Pe}=100$ and $\textrm{Pe}=30$ respectively, however due to the reduced conductance imposed by the hydrogel structures, the effective velocity for $\textrm{Pe}^{*}=4$ is $v=2$ $\mathrm{\mu m/s}$ and for $\textrm{Pe}^{*}=2$, $v=1$ $\mathrm{\mu m/s}$. In summary, it does not appear technically unfeasible to experimentally validate the existence of damped oscillations in active porous media.