1 Introduction
Advective–diffusive transport of passive or reactive scalars in confined environments (e.g. tubes and channels) is often accompanied by diffusive losses/gains through the confining walls. Examples of this phenomenon range from transport of nanoparticles in blood flow (Gentile, Ferrari & Decuzzi Reference Gentile, Ferrari and Decuzzi2008), growth of filamentous fungi (Heaton et al. Reference Heaton, López, Maini, Fricker and Jones2012a ,Reference Heaton, Obara, Grau, Jones, Nakagaki, Boddy and Fricker b ) and spread of contrast agents used in enhanced magnetic resonance imaging (Barnes, Quarles & Yankeelov Reference Barnes, Quarles and Yankeelov2014), through contaminant transport (Tang, Frind & Sudicky Reference Tang, Frind and Sudicky1981; Roubinet, de Dreuzy & Tartakovsky Reference Roubinet, de Dreuzy and Tartakovsky2012) and heat transfer (Martinez, Roubinet & Tartakovsky Reference Martinez, Roubinet and Tartakovsky2014; Gisladottir, Roubinet & Tartakovsky Reference Gisladottir, Roubinet and Tartakovsky2016) in fractured rocks, to microfluidics devices (Ling, Tartakovsky & Battiato Reference Ling, Tartakovsky and Battiato2016). We formulate our model, and interpret its results, in terms relevant to solute transport in tubes, but our analysis is directly transferable to other applications.
The relevant physicochemical phenomena are typically described by two sets of partial differential equations (PDEs) defined for a tube and the surrounding medium, which are coupled via continuity conditions at the tube’s wall. Numerical solutions of such problems are computationally demanding (often prohibitively so) due to discrepancy between the characteristic spatio-temporal scales of these PDEs (e.g. Gisladottir et al. Reference Gisladottir, Roubinet and Tartakovsky2016, and references therein). Previous analytical treatments of such problems relied on simplifying assumptions, such as steady-state flow and transport in a tube (Somers Reference Somers1912; Owen Reference Owen1925) or time-dependent transport in two-dimensional (Cartesian) channels of (semi)infinite extent (Tang et al. Reference Tang, Frind and Sudicky1981; Smith, Moore & Layton Reference Smith, Moore and Layton2003), which are not applicable for transient models of dispersion on networks.
We develop (semi-)analytical steady and transient one-dimensional solutions of advective–diffusive transport of solutes in a single tube with diffusive losses through the tube’s walls and first-order reaction rates in the surrounding medium. The time-dependent solutions are exact in the Laplace space, but rely on numerical inversion of Laplace transforms by means of the de Hoog algorithm (de Hoog, Knight & Stokes Reference de Hoog, Knight and Stokes1982; Hollenbeck Reference Hollenbeck1998). To facilitate the subsequent use of such solutions as building blocks in network models (e.g. Olufsen Reference Olufsen1999; Cousins & Gremaud Reference Cousins and Gremaud2012; Cousins, Gremaud & Tartakovsky Reference Cousins, Gremaud and Tartakovsky2013), we allow for the presence of general time-dependent boundary conditions at a tube’s inlet and outlet.
A mathematical model of solute transport in a tube and its invasion of, and degradation in, the surrounding medium is formulated in § 2. An analytical solution of this problem is presented in § 3 in terms of the Laplace-transformed solute concentrations in a tube (§ 3.1) and the surrounding medium (§ 3.2). These transforms are inverted analytically for both early times (§ 3.2.1) and steady state (§ 3.2.2). Results derived from our model are analysed in § 4. Finally, major conclusions are summarized in § 5.
2 Problem formulation
The analytical solutions derived below are applicable to a broad class of problems involving solute transport (described by a linear advection–diffusion–reaction equation) in a tube with diffusive losses through the tube’s walls.
2.1 Mathematical formulation
Consider a
$d$
-dimensional simulation domain
${\mathcal{D}}\subset \mathbb{R}^{d}$
that consists of two non-overlapping subdomains
${\mathcal{D}}_{1}$
and
${\mathcal{D}}_{2}$
such that
${\mathcal{D}}_{1}\cup {\mathcal{D}}_{2}={\mathcal{D}}$
and
${\mathcal{D}}_{1}\cap {\mathcal{D}}_{2}=\emptyset$
. Let
$\unicode[STIX]{x2202}{\mathcal{D}}$
and
$\unicode[STIX]{x1D6E4}$
, both in
$\mathbb{R}^{d-1}$
, denote the bounding surface of
${\mathcal{D}}$
and the interface separating
${\mathcal{D}}_{1}$
from
${\mathcal{D}}_{2}$
, respectively. The solute concentration
$C(\boldsymbol{x},t):{\mathcal{D}}\times \mathbb{R}^{+}\rightarrow \mathbb{R}^{+}$
is described by the solution of the system of linear advection–reaction–diffusion equations:
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}t}=-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{J}-\unicode[STIX]{x1D705}_{i}C,\quad \boldsymbol{J}=-D_{i}\unicode[STIX]{x1D735}C+\boldsymbol{u}_{i}C,\quad (\boldsymbol{x},t)\in {\mathcal{D}}_{i}\times (0,\infty ),\quad i=1,2.\end{eqnarray}$$
Here
$D_{i}>0$
are the diffusion coefficients,
$\boldsymbol{u}_{i}\in \mathbb{R}^{d}$
are the advective velocities,
$\boldsymbol{J}$
is the advection–diffusion flux, and
$\unicode[STIX]{x1D705}_{i}\geqslant 0$
are the reaction rate constants. These equations are coupled by enforcing continuity conditions at the interface
$\unicode[STIX]{x1D6E4}$
:
$$\begin{eqnarray}C(\boldsymbol{x}^{-},t)=C(\boldsymbol{x}^{+},t)\quad \text{and}\quad \boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{J}(\boldsymbol{x}^{-},t)=\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{J}(\boldsymbol{x}^{+},t),\quad \boldsymbol{x}_{\unicode[STIX]{x1D6E4}}\in \unicode[STIX]{x1D6E4}.\end{eqnarray}$$
Here
$\boldsymbol{n}(\boldsymbol{x}_{\unicode[STIX]{x1D6E4}})$
is the unit normal vector to
$\unicode[STIX]{x1D6E4}$
at point
$\boldsymbol{x}_{\unicode[STIX]{x1D6E4}}$
, and
$\boldsymbol{x}^{-}$
and
$\boldsymbol{x}^{+}$
indicate the limits of
$C$
and
$\boldsymbol{J}$
as
$\boldsymbol{x}\rightarrow \boldsymbol{x}_{\unicode[STIX]{x1D6E4}}$
from inside subdomains
${\mathcal{D}}_{1}$
and
${\mathcal{D}}_{2}$
, respectively. They are also subject to appropriate initial and boundary conditions on the domain boundary
$\unicode[STIX]{x2202}{\mathcal{D}}$
.
In the context of solute transport in a tube with radius
$a$
and length
$L$
(
$a\ll L$
), the tube is represented by domain
${\mathcal{D}}_{1}=\{(r,x):0\leqslant r<a,0<x<L\}$
and the surrounding medium by
${\mathcal{D}}_{2}=\{(r,x):a\leqslant r<\infty ,0<x<L\}$
(see figure 1). The advection–reaction–diffusion equations (2.1) take the form
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}t}=\frac{D_{1}}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(r\frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}r}\right)+D_{1}\frac{\unicode[STIX]{x2202}^{2}C}{\unicode[STIX]{x2202}x^{2}}-u(r)\frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}x}-\unicode[STIX]{x1D705}_{1}C,\quad (r,x)\in {\mathcal{D}}_{1},\end{eqnarray}$$
and
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}t}=\frac{D_{2}}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(r\frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}r}\right)+D_{2}\frac{\unicode[STIX]{x2202}^{2}C}{\unicode[STIX]{x2202}x^{2}}-\unicode[STIX]{x1D705}_{2}C,\quad (r,x)\in {\mathcal{D}}_{2}.\end{eqnarray}$$
The flow velocity
$u(r)$
is assumed to satisfy the Poiseuille (parabolic) law. The continuity conditions (2.2) at the interface
$\unicode[STIX]{x1D6E4}=\{(r,x):r=a,0\leqslant x\leqslant L\}$
become
$$\begin{eqnarray}C(a^{-},\cdot )=C(a^{+},\cdot ),\quad D_{1}\frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}r}(a^{-},\cdot )=D_{2}\frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}r}(a^{+},\cdot )\equiv -J_{m}.\end{eqnarray}$$
Two key quantities of interest are (i) total diffusive loss of solute through the tube’s walls and (ii) volume of the host medium affected by the solute.
Figure 1. Schematic representation of flow in a tube with diffusive losses to the surrounding medium.
2.2 Hydrodynamic dispersion approximation
Combined with the geometric constraint
$a\ll L$
, these quantities of interest justify approximating
$C(r,x,t)$
in the tube
${\mathcal{D}}_{1}$
with its cross-sectionally averaged counterpart,
$$\begin{eqnarray}C_{av}(x,t)=\frac{2}{a^{2}}\int _{0}^{a}C(r,x,t)r\,\text{d}r.\end{eqnarray}$$
It follows from (2.3a
) that the dynamics of
$C_{av}(x,\,t)$
satisfies approximately an advection–reaction–dispersion equation (Taylor Reference Taylor1953; Aris Reference Aris1956)
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}C_{av}}{\unicode[STIX]{x2202}t}=-\frac{2}{a}J_{m}+D\frac{\unicode[STIX]{x2202}^{2}C_{av}}{\unicode[STIX]{x2202}x^{2}}-v\frac{\unicode[STIX]{x2202}C_{av}}{\unicode[STIX]{x2202}x}-\unicode[STIX]{x1D705}_{1}C_{av},\quad 0<x<L,\end{eqnarray}$$
where
$v\equiv 2a^{-2}\!\int _{0}^{a}u(r)r\,\text{d}r$
is the average flow velocity,
$D=D_{1}+\unicode[STIX]{x1D6FC}v$
is the hydrodynamic dispersion coefficient with dispersivity
$\unicode[STIX]{x1D6FC}=a^{2}v/(48D_{1})$
, and
$J_{m}(z,t)$
is the (as yet unknown) diffusive flux from the tube into the surrounding medium. This equation is subject to the initial and boundary conditions
$$\begin{eqnarray}C_{av}(x,0)=c_{in}(x),\quad C_{av}(0,t)=c_{0}(t),\quad C_{av}(L,t)=c_{L}(t),\quad\end{eqnarray}$$
where
$c_{in}(x)$
is the (prescribed) initial concentration of solute in the tube, and
$c_{0}(t)$
and
$c_{L}(t)$
are the (prescribed) solute concentrations at the tube’s inlet (
$x=0$
) and outlet (
$x=L$
).
Interfacial continuity conditions (2.4) are replaced with a boundary condition for solute concentration in the medium, i.e. for reaction–diffusion equation (2.3b ),
$$\begin{eqnarray}C(a,x,t)=C_{av}(x,t).\end{eqnarray}$$
This is supplemented with the initial and boundary conditions
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle C(r,x,0)={\mathcal{C}}_{0},\quad C(\infty ,x,t)<\infty ,\\ \displaystyle \frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}x}(r,0,t)=0,\quad \frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}x}(r,L,t)=0.\end{array}\right\}\end{eqnarray}$$
2.3 Dimensionless formulation
Let us introduce dimensionless variables,
$$\begin{eqnarray}\hat{r}=\frac{r}{a},\quad \hat{x}=\frac{x}{L},\quad \hat{t}=\frac{tD_{2}}{a^{2}},\quad {\hat{C}}=\frac{C}{{\mathcal{C}}_{0}},\quad {\hat{C}}_{av}=\frac{C_{av}}{{\mathcal{C}}_{0}},\quad {\hat{c}}_{i}=\frac{c_{i}}{{\mathcal{C}}_{0}}\quad \text{for}~i=in,0,L,\quad\end{eqnarray}$$
and dimensionless parameters,
$$\begin{eqnarray}\unicode[STIX]{x1D700}=\frac{a}{L}\ll 1,\quad \unicode[STIX]{x1D6FC}_{D}=\frac{D}{D_{2}}>1,\quad Da_{1}=\frac{\unicode[STIX]{x1D705}_{1}L^{2}}{D},\quad Da_{2}=\frac{\unicode[STIX]{x1D705}_{2}a^{2}}{D_{2}},\quad Pe=\frac{vL}{D}.\end{eqnarray}$$
Then the boundary-value problem (BVP) for the medium, (2.3b
) and (2.8), takes a dimensionless form (
${\hat{c}}_{in,2}=1.0$
),
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}{\hat{C}}}{\unicode[STIX]{x2202}\hat{t}}=\frac{1}{\hat{r}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\hat{r}}\left(\hat{r}\frac{\unicode[STIX]{x2202}{\hat{C}}}{\unicode[STIX]{x2202}\hat{r}}\right)+\unicode[STIX]{x1D700}^{2}\frac{\unicode[STIX]{x2202}^{2}{\hat{C}}}{\unicode[STIX]{x2202}\hat{x}^{2}}-Da_{2}{\hat{C}},\quad 1<\hat{r}<\infty ,\quad 0<\hat{x}<1,\end{eqnarray}$$
$$\begin{eqnarray}{\hat{C}}(\hat{r},\hat{x},0)=1,\quad {\hat{C}}(1,\hat{x},\hat{t})={\hat{C}}_{av}(\hat{x},\hat{t}),\quad {\hat{C}}(\infty ,\hat{x},\hat{t})<\infty ,\end{eqnarray}$$
supplemented with boundary conditions
$\unicode[STIX]{x2202}{\hat{C}}/\unicode[STIX]{x2202}\hat{x}=0$
at
$\hat{x}=0$
and
$1$
. Likewise, a dimensionless form for the BVP in the tube, (2.6) and (2.7), is
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}{\hat{C}}_{av}}{\unicode[STIX]{x2202}\hat{t}}=\unicode[STIX]{x1D700}^{2}\unicode[STIX]{x1D6FC}_{D}\left(\frac{\unicode[STIX]{x2202}^{2}{\hat{C}}_{av}}{\unicode[STIX]{x2202}\hat{x}^{2}}-Pe\frac{\unicode[STIX]{x2202}{\hat{C}}_{av}}{\unicode[STIX]{x2202}\hat{x}}-Da_{1}{\hat{C}}_{av}\right)-\frac{2a}{D_{2}}J_{m},\quad 0<\hat{x}<1,\end{eqnarray}$$
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle {\hat{C}}_{av}(\hat{x},0)={\hat{c}}_{in}(\hat{x}),\quad {\hat{C}}_{av}(0,\hat{t})={\hat{c}}_{0}(\hat{t}),\\ \displaystyle {\hat{C}}_{av}(1,\hat{t})={\hat{c}}_{L}(\hat{t}).\end{array}\right\}\end{eqnarray}$$
In the remainder of this study we deal with the dimensionless quantities, but drop the hat
$\hat{~}$
to simplify the notation.
3 Analytical solutions
The BVPs (2.11) and (2.12) are solved below via Laplace transformation. The latter is defined as
$$\begin{eqnarray}\tilde{C}(\cdot ,s)=\int _{0}^{\infty }C(\cdot ,t)\text{e}^{-st}\,\text{d}t,\end{eqnarray}$$
where
$s\in \mathbb{C}$
is the Laplace variable.
3.1 Solute concentration in the medium
Our solution is relevant under conditions wherein the solute concentration gradients in the medium satisfy
$\unicode[STIX]{x2202}C/\unicode[STIX]{x2202}x<\unicode[STIX]{x2202}C/\unicode[STIX]{x2202}r$
. The leading-order (in
$\unicode[STIX]{x1D700}$
) approximation of (2.11) is
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}t}=\frac{1}{r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(r\frac{\unicode[STIX]{x2202}C}{\unicode[STIX]{x2202}r}\right)-Da_{2}C,\quad 1<r<\infty ,\end{eqnarray}$$
which is subject to the auxiliary conditions (2.11b ). The Laplace transform of the solution to the resulting BVP is given by (see appendix A)
$$\begin{eqnarray}\tilde{C}(r,x,s)=[\tilde{C}_{av}(x,s)-\unicode[STIX]{x1D6FD}^{-2}]\frac{\text{K}_{0}(\unicode[STIX]{x1D6FD}r)}{\text{K}_{0}(\unicode[STIX]{x1D6FD})}+\unicode[STIX]{x1D6FD}^{-2},\quad \unicode[STIX]{x1D6FD}(s)=\sqrt{s+Da_{2}},\end{eqnarray}$$
where
$\text{K}_{n}(\cdot )$
, with
$n=0,1,\ldots ,$
denotes the
$n$
th-order modified Bessel function.
3.2 Solute concentration in the tube
The Laplace transform,
$\tilde{J}_{m}$
, of the diffusive flux, from the tube into the surrounding medium, first defined in (2.6), is computed from (3.3) as
$$\begin{eqnarray}\tilde{J}_{m}(x,s)\equiv -D_{2}\frac{\unicode[STIX]{x2202}\tilde{C}}{\unicode[STIX]{x2202}r}(1,x,s)=D_{2}\unicode[STIX]{x1D6FD}[\tilde{C}_{av}(x,s)-\unicode[STIX]{x1D6FD}^{-2}]\frac{\text{K}_{1}(\unicode[STIX]{x1D6FD})}{\text{K}_{0}(\unicode[STIX]{x1D6FD})}.\end{eqnarray}$$
With this result, a Laplace-transformed solution of (2.12) is given by (appendix B)
$$\begin{eqnarray}\tilde{C}_{av}(x,s)=\unicode[STIX]{x1D6FE}+\frac{(\tilde{c}_{0}-\unicode[STIX]{x1D6FE})\text{e}^{\unicode[STIX]{x1D6FE}_{-}}-\tilde{c}_{L}+\unicode[STIX]{x1D6FE}}{\text{e}^{\unicode[STIX]{x1D6FE}_{-}}-\text{e}^{\unicode[STIX]{x1D6FE}_{+}}}\text{e}^{\unicode[STIX]{x1D6FE}_{+}x}+\frac{\tilde{c}_{L}-(\tilde{c}_{0}-\unicode[STIX]{x1D6FE})\text{e}^{\unicode[STIX]{x1D6FE}_{+}}-\unicode[STIX]{x1D6FE}}{\text{e}^{\unicode[STIX]{x1D6FE}_{-}}-\text{e}^{\unicode[STIX]{x1D6FE}_{+}}}\text{e}^{\unicode[STIX]{x1D6FE}_{-}x},\end{eqnarray}$$
where
$$\begin{eqnarray}\unicode[STIX]{x1D6FE}=\left[c_{in}+\frac{2}{\unicode[STIX]{x1D6FD}}\frac{\text{K}_{1}(\unicode[STIX]{x1D6FD})}{\text{K}_{0}(\unicode[STIX]{x1D6FD})}\right]\left[s+\unicode[STIX]{x1D700}^{2}\unicode[STIX]{x1D6FC}_{D}Da_{1}+2\unicode[STIX]{x1D6FD}\frac{\text{K}_{1}(\unicode[STIX]{x1D6FD})}{\text{K}_{0}(\unicode[STIX]{x1D6FD})}\right]^{-1}\end{eqnarray}$$
and
$$\begin{eqnarray}\unicode[STIX]{x1D6FE}_{\pm }=\frac{Pe}{2}\pm \frac{1}{2}\sqrt{Pe^{2}+\frac{4}{\unicode[STIX]{x1D700}^{2}\unicode[STIX]{x1D6FC}_{D}}\left[s+\unicode[STIX]{x1D700}^{2}\unicode[STIX]{x1D6FC}_{D}Da_{1}+2\unicode[STIX]{x1D6FD}\frac{\text{K}_{1}(\unicode[STIX]{x1D6FD})}{\text{K}_{0}(\unicode[STIX]{x1D6FD})}\right]}.\end{eqnarray}$$
The Laplace-transformed solution (3.5) is inverted numerically using the de Hoog algorithm (de Hoog et al. Reference de Hoog, Knight and Stokes1982; Hollenbeck Reference Hollenbeck1998). For early times and steady state, the inversion is carried out analytically in §§ 3.2.1 and 3.2.2.
3.2.1 Early-time solution
For small times
$t$
, i.e. large
$\text{Re}\,s$
or
$\text{Re}\,\unicode[STIX]{x1D6FD}$
, an asymptotic expansion of the modified Bessel functions is (Abramowitz & Stegun Reference Abramowitz and Stegun1984, equation (9.7.2))
$$\begin{eqnarray}\text{K}_{0}(\unicode[STIX]{x1D6FD})=\text{K}_{1}(\unicode[STIX]{x1D6FD})=\sqrt{\frac{\unicode[STIX]{x03C0}}{2\unicode[STIX]{x1D6FD}}}\,\text{e}^{-\unicode[STIX]{x1D6FD}}+O(\unicode[STIX]{x1D6FD}^{-3/2}\text{e}^{-\unicode[STIX]{x1D6FD}}),\end{eqnarray}$$
and (3.5b ) and (3.5c ) reduce to
$$\begin{eqnarray}\unicode[STIX]{x1D6FE}=\frac{c_{in}}{s}+O(s^{-3/2}),\quad \unicode[STIX]{x1D6FE}_{\pm }=\pm \frac{\sqrt{s}}{\unicode[STIX]{x1D700}\sqrt{\unicode[STIX]{x1D6FC}_{D}}}+O(s^{1/4}),\end{eqnarray}$$
so that
$\unicode[STIX]{x1D6FE}_{-}=-\unicode[STIX]{x1D6FE}_{+}$
. Hence,
$\exp (\unicode[STIX]{x1D6FE}_{-})\ll \exp (\unicode[STIX]{x1D6FE}_{+})$
and (3.5) is approximated by
$$\begin{eqnarray}\tilde{C}_{av}(x,s)=\unicode[STIX]{x1D6FE}+(\tilde{c}_{L}-\unicode[STIX]{x1D6FE})[\text{e}^{\unicode[STIX]{x1D6FE}_{+}(x-1)}-\text{e}^{-\unicode[STIX]{x1D6FE}_{+}(x+1)}]+(\tilde{c}_{0}-\unicode[STIX]{x1D6FE})[\text{e}^{-\unicode[STIX]{x1D6FE}_{+}x}-\text{e}^{\unicode[STIX]{x1D6FE}_{+}(x-2)}].\end{eqnarray}$$
For the constant boundary functions
$c_{0}$
and
$c_{L}$
, the analytical inversion of this Laplace-transformed solution rests on the relation
${\mathcal{L}}^{-1}\{\exp (a\sqrt{s})/s\}=1+\text{erf}[a/(2\sqrt{t})]$
, which leads to
$$\begin{eqnarray}\displaystyle C_{av}(x,t) & = & \displaystyle c_{in}+(c_{L}-c_{in})\left[\text{erf}\left(\frac{x-1}{2\unicode[STIX]{x1D700}\sqrt{\unicode[STIX]{x1D6FC}_{D}t}}\right)+\text{erf}\left(\frac{x+1}{2\unicode[STIX]{x1D700}\sqrt{\unicode[STIX]{x1D6FC}_{D}t}}\right)\right]\nonumber\\ \displaystyle & & \displaystyle -\,(c_{0}-c_{in})\left[\text{erf}\left(\frac{x}{2\unicode[STIX]{x1D700}\sqrt{\unicode[STIX]{x1D6FC}_{D}t}}\right)+\text{erf}\left(\frac{x-2}{2\unicode[STIX]{x1D700}\sqrt{\unicode[STIX]{x1D6FC}_{D}t}}\right)\right].\end{eqnarray}$$
3.2.2 Steady-state solution
The steady-state solution is obtained from (3.5) by invoking the final-value theorem,
$$\begin{eqnarray}C_{av}(x,\infty )=\mathop{\lim }_{s\rightarrow 0}s\tilde{C}_{av}(x,s).\end{eqnarray}$$
For constant boundary functions
$c_{0}$
and
$c_{L}$
, this yields
$$\begin{eqnarray}C_{av}(x,\infty )=\frac{c_{0}\text{e}^{\unicode[STIX]{x1D6FE}_{-}}-c_{L}}{\text{e}^{\unicode[STIX]{x1D6FE}_{-}}-\text{e}^{\unicode[STIX]{x1D6FE}_{+}}}\text{e}^{\unicode[STIX]{x1D6FE}_{+}x}+\frac{c_{L}-c_{0}\text{e}^{\unicode[STIX]{x1D6FE}_{+}}}{\text{e}^{\unicode[STIX]{x1D6FE}_{-}}-\text{e}^{\unicode[STIX]{x1D6FE}_{+}}}\text{e}^{\unicode[STIX]{x1D6FE}_{-}x},\end{eqnarray}$$
where
$$\begin{eqnarray}\unicode[STIX]{x1D6FE}_{\pm }=\frac{Pe}{2}\pm \sqrt{\frac{Pe^{2}}{4}+Da_{1}+2\frac{\sqrt{Da_{2}}}{\unicode[STIX]{x1D700}^{2}\unicode[STIX]{x1D6FC}_{D}}\frac{\text{K}_{1}(\sqrt{Da_{2}}\,)}{\text{K}_{0}(\sqrt{Da_{2}})}}.\end{eqnarray}$$
In addition to these solutions for a tube of finite length, appendix C contains a solution for a semi-infinite vessel. Such solutions can be used to analyse a single tube connected to an aggregated network with an unimpeded outflow condition.
4 Simulation results
For the sake of simplicity, in verifying our model we neglect the effects of the reaction terms, i.e. we set
$Da_{1}=Da_{2}=0$
. We define a dimensionless parameter
$\unicode[STIX]{x1D701}=a^{2}v/(LD_{2})$
, which is the ratio between the solute residence time in a tube,
$L/v$
, and the rate of diffusive losses through the tube’s wall,
$a^{2}/D_{2}$
. Consequently,
$\unicode[STIX]{x1D701}=Pe\,\unicode[STIX]{x1D700}^{2}$
. Figure 2 exhibits temporal snapshots of concentration profiles for
$\unicode[STIX]{x1D701}$
that spans two orders of magnitude,
$\unicode[STIX]{x1D701}=0.1$
and
$\unicode[STIX]{x1D701}=10.0$
. The dimensionless initial and boundary functions are set to
$c_{in}=0.1$
,
$c_{0}=1.0$
and, unless noted otherwise,
$c_{L}=0.1$
.
Figure 2. Temporal snapshots of dimensionless concentration profiles,
$C_{av}(x,\,\cdot )$
, computed with our analytical solution (3.5) (solid lines) and by solving numerically the full problem (2.11) and (2.12) (crosses), for (a)
$\unicode[STIX]{x1D701}=0.1$
and (b)
$\unicode[STIX]{x1D701}=10.0$
.
4.1 Model verification
The analytical solutions developed above represent the leading-order (in the radius-to-length ratio
$\unicode[STIX]{x1D700}$
) approximation of the full problem (see § 3.1), which neglects longitudinal variability of solute concentration in the surrounding medium. To verify the veracity of this approximation, we compare the inverse transform of
$\tilde{C}_{av}(x,s)$
with a numerical solution of the full two-dimensional BVPs (2.11) and (2.12). The software COMSOL was used to obtain the latter solution. Figure 2 demonstrates close agreement between the two solutions for both values of the dimensional parameter
$\unicode[STIX]{x1D701}$
. Similar agreement between the two solutions is observed for semi-infinite vessels as well (see appendix C).
For the prescribed (and constant) concentrations at the tube’s inlet and outlet, the solute concentrations increase with time. The rate of this increase is controlled by the dimensionless parameter
$\unicode[STIX]{x1D701}$
(see figure 2): the solute concentration increases faster for larger
$\unicode[STIX]{x1D701}$
, as the effects of solute advection in the tube increasingly dominate diffusive losses through the tube’s wall.
4.2 Flux through vessel wall
A key quantity of interest is the (dimensionless) total diffusive flux through a tube’s wall,
$$\begin{eqnarray}J_{tot}(t)=\frac{2a}{D_{2}c_{0}}\int _{0}^{L}J_{m}(x,t)\,\text{d}x.\end{eqnarray}$$
This is shown in figure 3 for a range of Péclet numbers. As Péclet number increases so does
${\hat{J}}_{tot}$
, but the ratio of solute diffused into the surrounding medium to solute supplied at the tube inlet decreases. This inverse relationship is due to the decreased residence time of solute in the tube. Large residence times allow greater diffusion through the tube wall. Furthermore, the maximum flux out of the tube occurs faster and with greater amplitude as
$Pe$
increases. That is because, as
$Pe$
increases, the rate at which the tube reaches its maximum solute concentration becomes progressively larger than the rate of solute consumption in the surrounding medium.
Figure 3. Total diffusive flux into the surrounding medium for several values of the Péclet number
$Pe$
with linear absorption kinetics.
As
$Pe$
increases, the
$J_{tot}(t)$
curves shift up and to the left. The leftward shift is largely due to our scaling of time
$t$
: increasing
$D_{2}$
, or decreasing
$a$
, causes the dimensionless time to grow, resulting in earlier peaks of
$J_{tot}(t)$
in figure 3. These peaks occur at the time at which the concentration front has reached the tube outlet, and the concentration gradient across the tube wall starts decaying to its equilibrium state. At the same time, higher values of
$D_{2}$
increase the diffusive flux through the tube wall, thereby increasing the amplitude of
$J_{tot}(t)$
.
Large Péclet numbers also result in sharper
$J_{tot}(t)$
profiles, while smooth
$J_{tot}(t)$
are characterized by low
$Pe$
. This is due to the shape of the concentration front transiting down the tube. Large-Péclet-number flows are dominated by the advective forces and the concentration front is steep. In low-Péclet-number flows the diffusive forces cause the concentration front to be more gradual. A steeper concentration front results in greater concentration gradients across the tube wall than observed with more gradual gradients.
Variation between the fluxes
$J_{tot}(t)$
computed with the finite and semi-infinite (for equal length) solutions, (3.5) and (C 2), is minimal. The finite-length-tube solution predicts that
$J_{tot}(t)$
is sustained for a longer period of time. That is because the imposed outflow condition causes more solute to exit through the vessel outlet and, hence, decreases the rate of contamination of the surrounding medium.
4.3 Effect of concentration gradient
Figure 4. Total diffusive flux into the surrounding medium for several values of the dimensionless outlet concentration
$c_{L}$
and for (a)
$Pe=5.0$
and (b)
$Pe=500.0$
, with linear adsorption kinetics.
The diffusive flux through the tube wall depends on the dimensionless concentration at the tube outlet. Figure 4 reveals that this dependence causes the
$J_{tot}(t)$
curve to shift up and to the left as
$c_{L}$
increases. The upward shift reflects the increased amount of solute available in the system. The leftward shift is due to the accelerated tube contamination as
$c_{L}$
increases, which is a byproduct of the first effect. The shape of the curves remains similar because the slope of the concentration front does not change significantly. This behaviour depends strongly on
$Pe$
because the advective flow diminishes the effects of the outflow boundary condition. As the Péclet number increases, the effects of the outlet boundary condition vanish; at
$Pe=500.0$
these effects are negligible.
It is worthwhile recalling that our dimensionless concentrations are defined relative to the inlet concentration
$c_{0}$
, such that
${\hat{c}}_{L}=c_{L}/c_{0}$
, with the hat dropped throughout much of the presentation for convenience.
4.4 Pulse boundary condition
A pulse boundary condition at the tube inlet is of relevance to a number of applications, including drug delivery, contaminant transport and heat transfer in liquid cooling or heating systems. We consider a finite pulse described by
$$\begin{eqnarray}\frac{c_{0}(t)}{{\mathcal{A}}_{0}}={\mathcal{H}}(t-T_{1})-{\mathcal{H}}(t-T_{2}),\end{eqnarray}$$
where
${\mathcal{A}}_{0}$
is the concentration of the pulse,
${\mathcal{H}}(\cdot )$
is the Heaviside function, and
$T_{1}$
and
$T_{2}$
are the times at which the pulse starts and ends. The Laplace transform of this boundary condition is
$$\begin{eqnarray}\tilde{c}_{0}(s)={\mathcal{A}}_{0}\frac{\exp (-T_{1}s)-\exp (-T_{2}s)}{s}.\end{eqnarray}$$
Figure 5. Temporal snapshots of dimensionless concentration profiles,
$C_{av}(x,\cdot )$
, for the pulse injection at the inlet (the concentration normalized with that of the pulse,
${\mathcal{A}}_{0}$
). All the dimensionless parameters are the same as above, except for
$c_{L}=c_{in}=0.1$
,
$T_{1}=0.0$
,
$T_{2}=0.01$
, and (a)
$\unicode[STIX]{x1D701}=10.0$
and (b)
$\unicode[STIX]{x1D701}=0.1$
. Results are presented as
$C_{av}/{\mathcal{A}}_{0}$
.
Figure 5 shows a solute pulse travelling down a tube. Recall that the dimensionless time
$t$
, defined in (2.9a–f
) and (2.10a–e
), rescales time with the dimensionless quantity
$\unicode[STIX]{x1D700}^{2}\unicode[STIX]{x1D6FC}_{D}$
. Increasing
$\unicode[STIX]{x1D700}^{2}\unicode[STIX]{x1D6FC}_{D}$
enhances the diffusive flux from the tube, causing the pulse to diminish faster. Increasing
$\unicode[STIX]{x1D700}^{2}\unicode[STIX]{x1D6FC}_{D}$
also decreases the width of the pulse and accelerates the movement of the pulse along the tube. The second phenomenon is a side effect of the first and is less prominent as the length of the tube decreases. It occurs when a solute diffuses back into the tube from the surrounding medium. The solute is absorbed by the surrounding medium at the front of the pulse; then the concentration gradient across the tube wall behind the pulse changes as the peak of the pulse passes the saturated medium, reversing the direction of the flux. This manifests itself in a pulse that travels slower and becomes wider. As
$\unicode[STIX]{x1D701}$
decreases, the pulse peak travels slower and experiences greater radial losses due to increased residence time.
5 Conclusions
We developed analytical solutions for advective–diffusive transport of chemically reactive solutes in a cylindrical tube with diffusive losses through its wall. The Laplace-transformed solutions are inverted analytically for early times and steady state, and numerically otherwise. Such solutions are relevant for a number of biological and engineered systems. Our analysis leads to the following major conclusions.
(i) Solute transport can be quantified by a dimensionless parameter
$\unicode[STIX]{x1D701}$
, which is the ratio of the solute’s residence time in a tube to the rate of diffusive losses through the tube wall.(ii) The impact of residence time in a tube on the rate of diffusive loss through the tube wall is ignored in network models, which assume a linear concentration distribution between a network’s nodes with a fixed solute loss through the tube walls (e.g. Heaton et al. Reference Heaton, López, Maini, Fricker and Jones2012a ,Reference Heaton, Obara, Grau, Jones, Nakagaki, Boddy and Fricker b ). The fixed-loss approximation might be reasonable for (sub)networks with geometrically similar tubes, but is inadequate for a heterogeneous network.
(iii) Our solutions for a finite-length tube satisfy Dirichlet boundary conditions with arbitrary time-dependent concentrations prescribed at the tube’s inlet and outlet. Consequently, they can be used as building blocks for models of solute transport in a bifurcating network in a manner similar to that used by Cousins & Gremaud (Reference Cousins and Gremaud2012) and Cousins et al. (Reference Cousins, Gremaud and Tartakovsky2013) to analyse network flows. In particular, in a given bifurcation, the outlet concentration of a source tube is prescribed as the inlet concentration of the respective downstream tubes. The advection–diffusion flux,
$\boldsymbol{J}$
, is also maintained across the bifurcation from source to downstream tubes. These two conservation conditions permit a solution to the concentration at the bifurcation as a function of the source inlet and downstream outlet conditions. Furthermore, these solutions can be linked for a network with a large number of tube generations, allowing a nodal network solution where only the source tube inlet and terminal tube outlet concentrations need to be known. However, it may be unrealistic to assume that the terminal tube outlet condition is knowable. In this case the use of the semi-infinite tube solution for the terminal tubes, which is provided in appendix C, obviates the need for the terminal tube outlet condition to be known.
In addition to analysing transport on networks, our solutions can be generalized in several ways. First, other state variables are consumed and/or generated inside and outside of a tube. Such scenarios can be accounted for by replacing the reactive term
$\unicode[STIX]{x1D705}C$
with its inhomogeneous counterpart
$\unicode[STIX]{x1D705}(C-C_{eq})$
, where
$C_{eq}$
is, for example, an equilibrium concentration.
Second, the infinite radial extent of the surrounding medium may be a low-fidelity representation of reality, which we used to model diffusive losses through the tube wall and subsequent contamination of the surrounding medium. In a setting considering numerous tubes embedded in a medium, a fair portion of solute ‘lost’ from the network may be reabsorbed by other tubes in the network. Our model can be modified to account for this effect by adding a sink term to (2.3b
) to represent a tube that may be absorbing solute or by replacing the semi-infinite domain
${\mathcal{D}}_{2}$
with a hollow cylinder whose outer surface is defined by an effective radius (average distance) of tube locations surrounding the tube of interest.
Appendix A. Laplace-transformed solution to BVP in the surrounding medium
The Laplace transformation of (3.2) yields an inhomogeneous Bessel equation (with
$\unicode[STIX]{x1D716}_{2}=\unicode[STIX]{x1D6FD}$
)
$$\begin{eqnarray}\frac{\text{d}^{2}\tilde{C}}{\text{d}r^{2}}+\frac{1}{r}\frac{\text{d}\tilde{C}}{\text{d}r}-\unicode[STIX]{x1D6FD}^{2}\tilde{C}=-1,\quad \unicode[STIX]{x1D6FD}(s)=\sqrt{s-Da_{2}}.\end{eqnarray}$$
Its general solution is the sum,
$$\begin{eqnarray}\tilde{C}(r;s)=C_{h}(r;s)+C_{p}(r;s),\end{eqnarray}$$
of a particular solution
$C_{p}$
to (A 1) and the general solution
$C_{h}$
to the corresponding homogeneous equation. The latter is given in terms of modified Bessel functions of the first,
$\text{I}_{0}(\cdot )$
, and second,
$\text{K}_{0}(\cdot )$
, kinds (Carslaw & Jaeger Reference Carslaw and Jaeger1959, p. 332):
$$\begin{eqnarray}C_{h}(r;s)=a_{1}\text{I}_{0}(\unicode[STIX]{x1D6FD}r)+a_{2}\text{K}_{0}(\unicode[STIX]{x1D6FD}r),\end{eqnarray}$$
where
$a_{1}$
and
$a_{2}$
are constants of integration. An inspection of (A 1) suggests its particular solution,
$C_{p}=1/\unicode[STIX]{x1D6FD}^{2}$
. Substituting these results into (A 2) yields
$$\begin{eqnarray}\tilde{C}(r,s)=a_{1}\text{I}_{0}(\unicode[STIX]{x1D6FD}r)+a_{2}\text{K}_{0}(\unicode[STIX]{x1D6FD}r)+\unicode[STIX]{x1D6FD}^{-2}.\end{eqnarray}$$
The constants
$a_{1}$
and
$a_{2}$
are obtained from the Laplace transforms of the boundary conditions (2.11b
),
$$\begin{eqnarray}\tilde{C}(1,s)=\tilde{C}_{av}(x,s),\quad \tilde{C}(\infty ,s)<\infty ,\end{eqnarray}$$
as
$a_{1}=0$
and
$a_{2}=(\tilde{C}_{av}-1/s)/\text{K}_{0}(\unicode[STIX]{x1D6FD})$
. The resulting Laplace-transformed solution for the tissue is given by (3.3).
Appendix B. Laplace-transformed solution to BVP in the tube
The Laplace transformation of (2.12), with the flux term
$\tilde{J}_{m}$
given by (3.4), yields
$$\begin{eqnarray}\frac{\text{d}^{2}\tilde{C}_{av}}{\text{d}x^{2}}-Pe\frac{\text{d}\tilde{C}_{av}}{\text{d}x}-F\tilde{C}_{av}=-G\end{eqnarray}$$
subject to
$$\begin{eqnarray}\tilde{C}_{av}(0,s)=\tilde{c}_{0}(s),\quad \tilde{C}_{av}(1,s)=\tilde{c}_{L}(s),\end{eqnarray}$$
where
$$\begin{eqnarray}F=\frac{1}{\unicode[STIX]{x1D700}^{2}\unicode[STIX]{x1D6FC}_{D}}\left[s-\unicode[STIX]{x1D700}^{2}\unicode[STIX]{x1D6FC}_{D}Da_{1}+2\unicode[STIX]{x1D6FD}\frac{\text{K}_{1}(\unicode[STIX]{x1D6FD})}{\text{K}_{0}(\unicode[STIX]{x1D6FD})}\right],\quad G=\frac{1}{\unicode[STIX]{x1D700}^{2}\unicode[STIX]{x1D6FC}_{D}}\left[c_{in}-\frac{2}{\unicode[STIX]{x1D6FD}}\frac{\text{K}_{1}(\unicode[STIX]{x1D6FD})}{\text{K}_{0}(\unicode[STIX]{x1D6FD})}\right].\end{eqnarray}$$
A solution of (B 1a
) is the sum of its particular solution,
$C_{av}^{(p)}$
, and the general solution,
$C_{av}^{(h)}$
, of the corresponding homogeneous equation,
$$\begin{eqnarray}\tilde{C}_{av}(x,s)=C_{av}^{(h)}(x,s)+C_{av}^{(p)}(x,s).\end{eqnarray}$$
Inspection of (B 1a
) reveals that a constant, e.g.
$\unicode[STIX]{x1D6FE}$
, is its particular solution, such that
$$\begin{eqnarray}C_{av}^{(p)}(x,s)\equiv \unicode[STIX]{x1D6FE}(s)=G/F.\end{eqnarray}$$
Since the characteristic equation of the homogeneous counterpart of (B 1a
), namely
$\unicode[STIX]{x1D706}^{2}-Pe\,\unicode[STIX]{x1D706}-F=0$
, has two real roots,
$$\begin{eqnarray}\unicode[STIX]{x1D6FE}_{\pm }=\frac{Pe}{2}\pm \frac{1}{2}\sqrt{Pe^{2}+4F},\end{eqnarray}$$
the homogeneous equation has the solution
$$\begin{eqnarray}C_{av}^{(h)}(x,s)=a_{1}\text{e}^{\unicode[STIX]{x1D6FE}_{+}x}+a_{2}\text{e}^{\unicode[STIX]{x1D6FE}_{-}x}.\end{eqnarray}$$
Substituting (B 4) and (B 6) into (B 3) yields
$$\begin{eqnarray}C_{av}(x,s)=a_{1}\text{e}^{\unicode[STIX]{x1D6FE}_{+}x}+a_{2}\text{e}^{\unicode[STIX]{x1D6FE}_{-}x}+\unicode[STIX]{x1D6FE}.\end{eqnarray}$$
It follows from (B 1b,c ) that
$$\begin{eqnarray}a_{1}=\frac{(\tilde{c}_{0}-\unicode[STIX]{x1D6FE})\text{e}^{\unicode[STIX]{x1D6FE}_{-}}-\tilde{c}_{L}+\unicode[STIX]{x1D6FE}}{\text{e}^{\unicode[STIX]{x1D6FE}_{-}}-\text{e}^{\unicode[STIX]{x1D6FE}_{+}}},\quad a_{2}=\frac{\tilde{c}_{L}-(\tilde{c}_{0}-\unicode[STIX]{x1D6FE})\text{e}^{\unicode[STIX]{x1D6FE}_{+}}-\unicode[STIX]{x1D6FE}}{\text{e}^{\unicode[STIX]{x1D6FE}_{-}}-\text{e}^{\unicode[STIX]{x1D6FE}_{+}}}.\end{eqnarray}$$
This gives rise to (3.5).
Appendix C. Solutions for a semi-infinite vessel
Consider (2.6) defined on the semi-infinite domain
$0\leqslant x\leqslant \infty$
and subject to the boundary conditions
$$\begin{eqnarray}C_{av}(0,t)=c_{0},\quad C_{av}(\infty ,t)=c_{in}.\end{eqnarray}$$
We use (3.3) to compute the Laplace transform of the (dimensional) diffusive solute flux from the tube into the surrounding medium. Following the steps described in appendices A and B, we obtain a Laplace-transformed solution
$$\begin{eqnarray}\tilde{C}_{av}(x,s)=\unicode[STIX]{x1D6FE}+(\tilde{c}_{0}-\unicode[STIX]{x1D6FE})\text{e}^{\unicode[STIX]{x1D6FE}_{-}x},\end{eqnarray}$$
where
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FE} & = & \displaystyle \left[c_{in}+\frac{2{\mathcal{C}}_{0}}{\unicode[STIX]{x1D6FD}\sqrt{D_{2}}}\frac{\text{K}_{1}(\unicode[STIX]{x1D6FD})}{\text{K}_{0}(\unicode[STIX]{x1D6FD})}\right]\left[s-\unicode[STIX]{x1D705}_{1}+\frac{2\sqrt{D_{2}}}{a^{2}}\frac{\text{K}_{1}(\unicode[STIX]{x1D6FD})}{\text{K}_{0}(\unicode[STIX]{x1D6FD})}\right]^{-1},\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FE}_{-} & = & \displaystyle \frac{v}{2D}-\frac{1}{2}\sqrt{\frac{v^{2}}{D^{2}}+\frac{4}{D}\left[s-\unicode[STIX]{x1D705}_{1}+\frac{2}{a^{2}}\frac{\text{K}_{1}\left(\unicode[STIX]{x1D6FD}\right)}{\text{K}_{0}(\unicode[STIX]{x1D6FD})}\right]},\end{eqnarray}$$
and
$\unicode[STIX]{x1D6FD}=a\sqrt{(s-\unicode[STIX]{x1D705}_{2})/D_{2}}$
.
The early-time solution is obtained analytically, for constant
$c_{0}$
, by using (3.6) to obtain an asymptotic expansion of (C 2) and inverting the resulting expression:
$$\begin{eqnarray}C_{av}(x,t)=c_{0}+c_{in}\,\text{erf}\left(\frac{x}{2\unicode[STIX]{x1D700}\sqrt{\unicode[STIX]{x1D6FC}_{D}t}}\right).\end{eqnarray}$$
The corresponding steady-state concentration is obtained from (C 2) by applying the finite-value theorem,
$$\begin{eqnarray}C_{av}(x,\infty )=c_{0}\text{e}^{\unicode[STIX]{x1D6FE}_{-}x}.\end{eqnarray}$$
Numerical inversion of (C 2),
$C_{av}(x,t)$
, is compared with a numerical solution of the full two-dimensional boundary-value problems (2.11) and (2.12) in figure 6. The close agreement between the two solutions demonstrates the adequacy of the assumption leading to (3.2) for the full range of interesting model parameters. This agreement improves as
$\unicode[STIX]{x1D701}$
increases.
