Steady state diffusion in tubular structures: Assessment of one-dimensional models

Abstract

Steady-state diffusion in long axisymmetric structures is considered. The goal is to assess one-dimensional approximations by comparing them with axisymmetric eigenfunction expansions. Two problems are considered in detail: a finite tube with one end that is partly absorbing and partly reflecting; and two finite coaxial tubes with different cross-sectional radii joined together abruptly. Both problems may be modelled using effective boundary conditions, containing a parameter known as the trapping rate. We show that trapping rates depend on the lengths of the finite tubes (and that they decay slowly as these lengths increase) and we show how trapping rates are related to blockage coefficients, which are well known in the context of potential flow along tubes of infinite length.

1 Introduction

Problems involving wave propagation or diffusion in long rigid tubes of slowly varying cross-section have an extensive literature, with associated approximations. One thinks of the Webster horn equation or the Fick–Jacobs equation, which are one-dimensional partial differential equations, in which the unknown function depends on the longitudinal coordinate z and time t; in steady situations (no dependence on t), ordinary differential equations are obtained. These approximations are not expected to be useful when the geometry changes abruptly.

What can we do if the slowly varying assumption is violated? This happens if there is an abrupt change in the tube’s diameter at $z=0$ or if the tube contains an object such as a thin disc in the plane $z=0$ . (More complicated situations are readily imagined, but the two examples mentioned are sufficient for our purposes.) One possibility is to connect the solution on the two sides of $z=0$ using certain interface conditions, while still retaining a one-dimensional model. This is an old idea, in which the effects of a boundary with various structures (such as perforations or decorations) are smeared out (homogenised) and replaced by an effective boundary or interface condition. These conditions involve parameters that are supposed to encapsulate properties of the boundary structures; these parameters have to be calculated or estimated by a separate procedure. In the diffusion literature, the approach just described is known as boundary homogenisation. Quoting from [4]: ‘The key step of our analysis is the replacement of the initial three-dimensional problem by an effective one-dimensional one. This replacement is not rigorous. Therefore, our analytical results do not provide an exact solution to the problem.’ This provides one motivation for our work.

Another possibility is to tackle the fully three-dimensional problem using eigenfunction expansions. Thus, if we assume that the tube has a circular cross-section, we can use cylindrical polar coordinates $(r, \theta, z)$ and the method of separation of variables. We shall use this method for two specific problems and then compare with one-dimensional approximations.

We shall give detailed consideration of two steady-state axisymmetric problems involving Laplace’s equation, $\nabla^2\phi=0$ , and tubular geometries. In the first problem (the ‘finite-tube problem’, Section 2), the rigid tube has finite length ( $0<z<\ell$ ) with $\phi=\phi_0$ (a constant) at one end ( $z=\ell$ ) and mixed conditions at the other ( $z=0$ , where $\partial {\phi}/\partial {z}=0$ for $0\leq r<b$ and $\phi=0$ for $b<r<a$ ). When boundary homogenisation is used, the goal is to calculate $\psi(z)\simeq\phi$ with $\psi''(z)=0$ . To be more precise, we may define

(1.1) \begin{equation} \psi(z)=\frac{1}{|{\mathcal A}(z)|}\int_{{\mathcal A}(z)}\phi(r,\theta,z)r\,{\rm d} r\,{\rm d}\theta,\end{equation}

where ${\mathcal A}(z)$ is the cross-section at z and $|{\mathcal A}(z)|$ is the area of that cross-section. The mixed boundary conditions at $z=0$ are replaced by a single (Robin) boundary condition,

(1.2) \begin{equation} \psi'(0)=\kappa \psi(0),\end{equation}

where $\kappa$ is a constant. In the diffusion literature, $\kappa$ is known as a trapping rate. Particles diffusing in the tube are partially absorbed at $z=0$ and partially reflected; $\kappa$ characterises the absorption rate.

Determining $\kappa$ is a separate and non-trivial task. For example, if we switch the mixed conditions at $z=0$ , so that $\phi=0$ for $0\leq r<b$ (representing an absorbing disc) and $\partial {\phi}/\partial {z}=0$ for $b<r<a$ (a reflecting annulus), there is the approximation given in [2], obtained there using Brownian dynamics simulations. Such approximations do not depend on $\ell$ , the length of the tube. A second motivation for our work was to quantify the dependence of $\kappa(\ell)$ on $\ell$ . We shall see that

(1.3) \begin{equation} \kappa(\ell)=\kappa_\infty +O\left(\ell^{-1}\right)\quad\mbox{as $\ell\to\infty$,}\end{equation}

where $\kappa_\infty$ comes from solving a potential flow problem for a tube of infinite length. In fact, we shall also see that $\kappa_\infty=C^{-1}$ , where C is the blockage coefficient for the potential flow. The slow decay with $\ell$ in (1.3) is surprising because one might expect to see exponential decay. We also obtain an estimate of the $O(\ell^{-1})$ term; it turns out to be simply $-\ell^{-1}$ (see (2.26)), a result obtained by applying an iterative method to an associated integral equation.

Much is known about blockage coefficients, such as how to compute them and how to estimate them for obstacles of various shapes in an infinite tube. For details and references, see [14, 15]. Connections between C and $\kappa$ are expected to find further applications.

In the second problem (the ‘junction problem’, Section 3), two coaxial rigid tubes are joined at $z=0$ . The tubes have different radii and different finite lengths. The join is effected using a rigid annulus, and there are Dirichlet boundary conditions at the two ends of the structure ( $\phi=\phi_0$ at one end and $\phi=0$ at the other). This problem was studied by Berezhkovskii et al. [1] (using boundary homogenisation, see Section 3.1) and by Kalinay and Percus [10] (using eigenfunction expansions). The junction problem is more complicated than the finite-tube problem: instead of mixed conditions at one end of the tube, we now have a mixture of boundary and continuity conditions at $z=0$ . It turns out that (1.2) is to be replaced by a condition requiring that $\psi'(0)$ is proportional to the jump in $\psi$ across $z=0$ ; the constant of proportionality is essentially $\kappa$ (see (3.8)), and it is found to decay slowly, similar to (1.3); see (3.32) for more details. We also relate $\kappa$ to the blockage coefficient for potential flow along two coaxial semi-infinite tubes, joined at $z=0$ , a problem that is studied in Section 3.2.

There are some concluding remarks in Section 4.

2 Finite-tube problem

Consider the following axisymmetric problem for a rigid tube of finite length $\ell$ . The potential $\phi(r,z)$ satisfies $\nabla^2\phi=0$ for $0<z<\ell$ , $0\leq r<a$ together with the following conditions:

(2.1a) \begin{align}& \phi(r,\ell)=\phi_0,\quad 0\leq r<a,\end{align}

(2.1b) \begin{align}& \partial {\phi}/\partial {z}=0,\quad z=0,\quad 0\leq r<b,\end{align}

(2.1c) \begin{align}& \phi(r,0)=0,\quad b<r<a,\end{align}

(2.1d) \begin{align}&\partial {\phi}/\partial {r}=0,\quad r=a, \quad 0<z<\ell.\end{align}

Here, $\phi_0$ is a constant. (Evidently, the problem could be scaled so that $\phi_0=1$ and $a=1$ , if desired.) A complementary problem (with the mixed conditions at $z=0$ interchanged, so that $\phi=0$ on the disc of radius b and $\partial {\phi}/\partial {z}=0$ on the annulus $b<r<a$ ) is studied in [4, 5].

The main quantity of interest is the flux

(2.2) \begin{equation} {\mathcal J}=2\pi \int_0^a \frac{\partial {\phi}}{\partial {z}}\,r\,{\rm d} r.\end{equation}

Green’s theorem shows that ${\mathcal J}$ does not depend on z, with $0<z<\ell$ , but it does depend on $\ell$ : we write ${\mathcal J}(\ell)$ . In particular, using (2.1b),

(2.3) \begin{equation} {\mathcal J}(\ell)= 2\pi \int_b^a \left.\frac{\partial {\phi}}{\partial {z}}\right|_{z=0}r\,{\rm d} r,\end{equation}

an integral of $\partial {\phi}/\partial {z}$ over the annular region at $z=0$ . Alternatively, an application of Green’s theorem in the finite tube to $\phi$ and $z-\ell$ gives

(2.4) \begin{equation} {\mathcal J}(\ell)=\pi a^2 \frac{\phi_0}{\ell} -\frac{2\pi}{\ell} \int_0^b \phi(r,0)r\,{\rm d} r,\end{equation}

giving ${\mathcal J}(\ell)$ in terms of an integral of $\phi$ over the disc at $z=0$ .

2.1 Boundary homogenisation

Suppose we approximate $\phi(r,z)$ by a function of z only, $\psi(z)$ ; see (1.1). Laplace’s equation becomes $\psi''(z)=0$ , so that $\psi(z)=Az+B$ , with arbitrary constants A and B. The mixed boundary conditions at $z=0$ are replaced by (1.2), $\psi'(0)=\kappa\psi(0)$ , where $\kappa$ is a constant. Combining this condition with $\psi(\ell)=\phi_0$ , we obtain

(2.5) \begin{equation} \psi(z)=\phi_0\,\frac{1+\kappa z}{1+\kappa\ell} \quad\mbox{and}\quad {\mathcal J}(\ell)\simeq \pi a^2 \psi' =\frac{\pi a^2\kappa\phi_0}{1+\kappa\ell}.\end{equation}

This calculation is very simple but, of course, we do not know $\kappa$ : the problem of calculating $\phi$ has been replaced by the problem of calculating $\kappa$ .

2.2 Rigid disc in an infinite tube

Returning to problem (2.1), suppose we extend $\phi(r,z)$ into $-\ell<z<0$ as an odd function of z, $\phi(r,-z)=-\phi(r,z)$ , $0<z<\ell$ . Then we see that the problem is equivalent to some kind of potential flow past a thin rigid disc of radius b in a tube of cross-sectional radius a. Note that extending $\phi$ in this way means that (2.1c) is satisfied automatically.

The simplest problem of this kind comes by letting $\ell\to\infty$ , with uniform flow along an infinite tube. For this problem, the conditions $\phi(r,\pm\ell)=\pm\phi_0$ (see (2.1a)) are replaced by

(2.6) \begin{equation} \phi(r,z)= U(z\pm C)+o(1)\quad\mbox{as $z\to\pm\infty$,}\end{equation}

where U is the constant speed of the flow far from the disc and C is the blockage coefficient. (The o(1) terms in (2.6) are exponentially small.) We shall see later that $\kappa$ in (1.2) and C in (2.6) are related.

The problem of uniform flow past a rigid disc has been studied previously [18, 17, 14, 15]. There are approximations when $b/a$ is small and when the gap between the disc and the tube is small ( $b/a\simeq 1$ ) [17, 14]. Here, we begin with a standard method (that could be used to compute C numerically) but then we introduce approximations that lead to analytical estimates.

Two exact formulas are worth noting. Hurley’s formula [9], [14, equation (31)] gives

(2.7) \begin{equation} C=\frac{2}{a^2U}\int_0^b \phi(r,0+)\,r\,{\rm d} r.\end{equation}

We also have flux conservation through the gap,

(2.8) \begin{equation} 2\pi \int_b^a \frac{\partial {\phi}}{\partial {z}}(r,0+)\,r\,{\rm d} r=\pi a^2 U.\end{equation}

Now, as it is sufficient to consider $z>0$ , we write

(2.9) \begin{equation} \phi(r,z)=U(z+C) +Ua\sum_{n=1}^\infty c_n J_0(\lambda_nr){\rm e}^{-\lambda_nz},\quad 0\leq r<a,\ z>0,\end{equation}

where $\lambda_n$ are positive solutions of $J_1(\lambda_na)=0$ and $J_n$ are Bessel functions. The (dimensionless) coefficients $c_n$ and $C/a$ are to be found. The representation (2.9) ensures that $\nabla^2\phi=0$ and that (2.1d) and (2.6) are satisfied. Applying (2.1b) and (2.1c) gives

(2.10a) \begin{align}& \sum_{n=1}^\infty c_n (\lambda_na) J_0(\lambda_nr)=1, \quad 0\leq r<b, \end{align}

(2.10b) \begin{align} & \frac{C}{a}+\sum_{n=1}^\infty c_n J_0(\lambda_nr)=0, \quad b<r<a.\end{align}

These are dual series equations for $c_n$ and C [19], [6, Example 3.3.1]. For a numerical treatment, see [17]. Also, if we multiply (2.10b) by r and then integrate from $r=b$ to $r=a$ , we obtain

(2.11) \begin{equation} \frac{C}{a}=\frac{2ab}{a^2-b^2}\sum_{n=1}^\infty c_n \frac{J_1(\lambda_nb)}{\lambda_na}\end{equation}

after use of

(2.12) \begin{equation} \int_p^q J_0(\lambda_nr)\,r\,{\rm d} r =\frac{q}{\lambda_n}\,J_1(\lambda_nq) -\frac{p}{\lambda_n}\,J_1(\lambda_np).\end{equation}

The same formula for C, (2.11), is obtained by substituting (2.9) in (2.7).

Instead of trying to solve (2.10), let $V(r)=\partial {\phi}/\partial {z}$ at $z=0$ for $b<r<a$ ; this is the flow speed in the gap. Thus, using (2.1b) and (2.9),

\begin{equation*} U-U\sum_{n=1}^\infty c_n (\lambda_n a)J_0(\lambda_nr)= \left\{\begin{array}{cl} 0, & 0\leq r<b,\\ \\[-7pt] V(r), & b<r<a. \end{array}\right.\end{equation*}

This is a Dini–Bessel series. Orthogonality gives

(2.13) \begin{equation} \pi a^2 U=2\pi\int_b^a V(r)r\,{\rm d} r\end{equation}

(which is (2.8)) and

(2.14) \begin{equation} c_n=-\frac{2}{U\lambda_na^3 J_0^2(\lambda_na)}\int_b^a V(r)J_0(\lambda_nr)r\,{\rm d} r,\end{equation}

using [16, 10.22.4 and 10.22.5]

(2.15) \begin{equation} \int_0^a J_0(\lambda_m r)J_0(\lambda_nr)\,r\,{\rm d} r= \frac{a^2}{2}\,J_0^2(\lambda_na)\,\delta_{mn}.\end{equation}

Substituting for $c_n$ from (2.14) in (2.10a) leads to an integral equation for V(r). This approach will be pursued later when we return to the finite-tube problem in Section 2.3. However, for the moment, we note that our calculations are exact, but we do not know V.

Let us make a plausible approximation for V. The simplest possibility is to take $V(r)=U_0$ , a constant; from (2.13), $U_0=Ua^2/(a^2-b^2)$ . Then (2.14) gives the approximation

\begin{equation*} c_n =\frac{2bU_0 J_1(\lambda_nb)}{U\lambda_n^2a^3 (a^2-b^2)\,J_0^2(\lambda_na)}\end{equation*}

which, when substituted in (2.11) gives the approximation

(2.16) \begin{equation} C\simeq C_{{\rm disc}}^{{\rm app}}=\frac{4a^3b^2}{(a^2-b^2)^2} \,{\mathbb S}_0(b/a)\end{equation}

where

(2.17) \begin{equation} {\mathbb S}_0(X)= \sum_{n=1}^\infty \frac{J_1^2(\lambda_naX)}{(\lambda_na)^3\,J_0^2(\lambda_na)} =\sum_{n=1}^\infty \frac{J_1^2(j_nX)}{j_n^3J_0^2(j_n)}\end{equation}

and $j_n=\lambda_na$ (so that $J_1(j_n)=0$ ). (The notation ‘app’ in (2.16) indicates ‘approximation’.) The terms in the series ${\mathbb S}_0$ decay as $n^{-3}$ ; as far as we know, ${\mathbb S}_0$ cannot be summed in closed form (although it can be written as an integral [13, equation (16)]). It will appear again later; see (2.24). A generalisation, ${\mathbb S}_q$ will also appear later; see (3.14) and Appendix A.

Although the approximation $V(r)\simeq U_0$ may seem crude, it can be refined, taking into account the fact that V(r) is infinite at the sharp edge of the disc ( $r=b$ ). For details, with an application to a disc that almost fills the tube (small gap, $b/a\simeq 1$ ), see [15].

2.3 Finite-tube problem: Estimating $\kappa$

2.3.1 An integral equation

Let us return to the finite-tube problem (2.1). Look for a solution in the form

(2.18) \begin{equation} \frac{\phi(r,z)}{\phi_0} =1+c_0\left(\frac{z}{\ell}-1\right) +\sum_{n=1}^\infty c_n J_0(\lambda_nr)\sinh{\{\lambda_n(z-\ell)\}}, \quad 0\leq r<a,\ 0<z<\ell.\end{equation}

This representation satisfies (2.1a) and (2.1d). From (2.2), the flux is given by

(2.19) \begin{equation} {\mathcal J}(\ell)=\pi a^2\phi_0 c_0/\ell.\end{equation}

Applying the boundary conditions at $z=0$ , (2.1b) and (2.1c), we obtain dual series equations,

(2.20a) \begin{align}c_0 +\sum_{n=1}^\infty c_n J_0(\lambda_nr) (\lambda_n\ell) \cosh{\lambda_n\ell} &=0,\quad 0\leq r<b, \end{align}

(2.20b) \begin{align}c_0 +\sum_{n=1}^\infty c_n J_0(\lambda_nr) \sinh{\lambda_n\ell} &=1,\quad b<r<a.\end{align}

Proceeding as in Section 2.2, let $V(r)=\partial {\phi}/\partial {z}$ at $z=0$ for $b<r<a$ . Then, using (2.1b) and (2.18), we obtain

\begin{equation*} c_0 +\sum_{n=1}^\infty c_n J_0(\lambda_nr)(\lambda_n\ell)\cosh{\lambda_n\ell} =\left\{\begin{array}{c@{\quad}l} 0, & 0\leq r<b,\\ \\[-7pt] (\ell/\phi_0) V(r), & b<r<a. \end{array}\right.\end{equation*}

Orthogonality gives

(2.21) \begin{equation} c_0=\frac{2\ell}{\phi_0a^2}\int_b^a V(r)r\,{\rm d} r =\frac{\ell}{\pi a^2\phi_0}{\mathcal J}(\ell)\end{equation}

and

\begin{equation*}c_n \lambda_n\ell \cosh{\lambda_n\ell}=\frac{2\ell}{\phi_0a^2J_0^2(\lambda_na) }\int_b^a V(r)J_0(\lambda_nr)r\,{\rm d} r\end{equation*}

using (2.15). Substitution in (2.20b) gives an integral equation for V,

(2.22) \begin{equation} 1= \frac{2\ell}{\phi_0a^2}\int_b^a V(s)\bigg\{ 1 +\frac{a}{\ell}\sum_{n=1}^\infty \frac{J_0(\lambda_nr)J_0(\lambda_ns)} {(\lambda_na) \,J_0^2(\lambda_na)} \tanh{\lambda_n\ell} \bigg\} s\,{\rm d} s,\quad b<r<a.\end{equation}

This reduction to an integral equation is exact. However, we shall approximate because we are only interested in calculating $c_0$ (see (2.21)) and, moreover, we shall assume that $\ell/a\gg 1$ .

2.3.2 Approximation for large $\ell$

For large $\ell$ , $\tanh{\lambda_n\ell}\simeq 1$ , so we can discard the sum term in (2.22) (because of the factor $a/\ell$ ). Thus, at leading order,

(2.23) \begin{equation} 1\simeq \frac{2\ell}{\phi_0a^2}\int_b^a V(s)s\,{\rm d} s=c_0,\end{equation}

which gives ${\mathcal J}(\ell)\simeq \pi a^2\phi_0/\ell$ . This very simple approximation does not depend on b, although it is exact when $b=0$ .

In order to improve our approximations, we use an iterative approach. Thus, suppose that $V(r)=V_0/\ell+V_1/\ell^2+\cdots$ , where $V_0$ is a constant defined by $\phi_0a^2=2V_0\int_b^ar\,{\rm d} r$ , which gives $V_0=\phi_0a^2/(a^2-b^2)$ . Substitution in (2.22) then gives

\begin{equation*} 1\simeq \frac{2\ell}{\phi_0a^2}\int_b^a \left(\frac{V_0}{\ell}+\frac{V_1}{\ell^2}\right)s\,{\rm d} s +\frac{2\ell}{\phi_0a^2}\int_b^a \frac{V_0a}{\ell^2} \sum_{n=1}^\infty \frac{J_0(\lambda_nr)J_0(\lambda_ns)} {(\lambda_na) \,J_0^2(\lambda_na)} s\,{\rm d} s,\quad b<r<a.\end{equation*}

Rearranging

\begin{align*} c_0&=1 -\frac{2V_0}{\phi_0a\ell} \sum_{n=1}^\infty \frac{J_0(\lambda_nr)} {(\lambda_na) \,J_0^2(\lambda_na)} \int_b^a J_0(\lambda_ns)s\,{\rm d} s\\[3pt] &=1 +\frac{2ba^2}{\ell(a^2-b^2)} \sum_{n=1}^\infty \frac{J_0(\lambda_nr)J_1(\lambda_nb)} {(\lambda_na)^2 \,J_0^2(\lambda_na)}.\end{align*}

But the right-hand side is a function of r, so average over the annulus: multiply by $2r/(a^2-b^2)$ and integrate between $r=b$ and $r=a$ :

(2.24) \begin{align}c_0&=1 +\frac{4ba^2}{\ell(a^2-b^2)^2} \sum_{n=1}^\infty\frac{J_1(\lambda_nb)} {(\lambda_na)^2 \,J_0^2(\lambda_na)}\int_b^aJ_0(\lambda_nr)r\,{\rm d} r\nonumber\\[3pt] &=1 -\frac{4b^2a^3}{\ell(a^2-b^2)^2} \sum_{n=1}^\infty\frac{J_1^2(\lambda_nb)} {(\lambda_na)^3 \,J_0^2(\lambda_na)}.\end{align}

The series seen here is ${\mathbb S}_0(b/a)$ , defined by (2.17). It occurred in our approximation to the blockage coefficient C for potential flow past a thin rigid disc in a tube of infinite length, $C_{{\rm disc}}^{{\rm app}}$ , given by (2.16). If we accept that approximation, we obtain the relation

(2.25) \begin{equation} c_0=1-C_{{\rm disc}}^{{\rm app}}/\ell\end{equation}

and then (2.19) gives an approximation for the flux ${\mathcal J}(\ell)$ .

2.3.3 Estimating $\kappa$

We have two estimates for ${\mathcal J}$ , namely

\begin{equation*} {\mathcal J}(\ell)=\frac{\pi a^2\phi_0}{\ell}\left(1-\frac{C_{{\rm disc}}^{{\rm app}}}{\ell}\right) \quad\mbox{and}\quad {\mathcal J}(\ell) =\frac{\pi a^2\phi_0\kappa}{1+\kappa\ell};\end{equation*}

see (2.5) for the second. Equating these and solving for $\kappa$ , we obtain

(2.26) \begin{equation} \kappa(\ell)=\frac{1}{C_{{\rm disc}}^{{\rm app}}} -\frac{1}{\ell}.\end{equation}

Recall that $C_{{\rm disc}}^{{\rm app}}$ pertains to an infinite tube whereas $\kappa(\ell)$ is for a finite tube of length $\ell$ . We see that the difference decays slowly as $\ell\to\infty$ .

The observation that $\kappa(\infty)=C^{-1}$ is to be expected. If we approximate $\phi(r,z)$ for an infinite pipe (with far field given by (2.6)) by a one-dimensional function $\psi(z)$ , we must have $\psi(z)=U(z+C)$ . Suppose the rigid object in the tube is confined to the region $|z|<h$ for some h, and then apply the condition $\psi'(z_0)=\kappa_\infty^0\psi(z_0)$ (compare with (1.2)), where $z_0\geq h$ ; we have written $\kappa_\infty^0$ to indicate that we are representing the solution in an infinite pipe with a boundary condition applied at $z_0$ . Doing this gives $\kappa_\infty^0=1/(C+z_0)$ . (For a thin disc, we can take $z_0=h=0$ .) This simple calculation shows that we should expect $\kappa$ to depend on the location at which the boundary condition $\psi'=\kappa\psi$ is imposed. Indeed, if we use a different location, using $\psi'(z_1)=\kappa_\infty^1\psi(z_1)$ , we obtain $\kappa_\infty^1=1/(C+z_1)$ and then, eliminating C, we obtain

(2.27) \begin{equation} \frac{1}{\kappa_\infty^0}-\frac{1}{\kappa_\infty^1}=z_0-z_1.\end{equation}

We conclude that, unlike the blockage coefficient C, $\kappa$ is not dependent solely on the shape of the rigid object in the tube.

2.4 Finite-tube problem: The one-dimensional approximation $\psi$

Let us construct $\psi$ from $\phi$ , using (1.1), which reduces to

\begin{equation*} \psi(z)=\frac{2}{a^2}\int_0^a \phi(r,z)r\,{\rm d} r=\phi_0\left\{1+c_0\left(\frac{z}{\ell}-1\right) \right\},\end{equation*}

an exact formula, after using (2.18). Evidently, $\psi(\ell)=\phi_0$ whereas imposition of $\psi'(0)=\kappa\psi(0)$ , (1.2), gives a relation between $\kappa$ and $c_0$ ,

(2.28) \begin{equation} \kappa=\frac{c_0}{\ell(1-c_0)}.\end{equation}

Now, the leading-order estimate for $c_0$ , (2.23), is $c_0\simeq 1$ , and so the exact formula (2.28) does not provide a useful estimate of $\kappa$ . However, if we use the refined estimate for $c_0$ , (2.25), in (2.28), we recover the formula (2.26).

3 Junction problem

Two finite rigid coaxial tubes join abruptly at $z=0$ . The wider tube has length $\ell_{{w}}$ and cross-sectional radius a. The narrower tube has length $\ell_{{n}}$ and cross-sectional radius $b\leq a$ . The potential $\phi$ satisfies $\partial {\phi}/\partial {r}=0$ on the tube walls together with $\partial {\phi}/\partial {z}=0$ on the annular region, $b<r<a$ , $z=0-$ . Dirichlet boundary conditions are imposed at the two ends of the structure. Thus, we have $\nabla^2\phi=0$ inside the structure, together with the following boundary conditions:

(3.1a) \begin{align}& \phi(r,\ell_{{n}})=\phi_0,\quad 0\leq r<b, \end{align}

(3.1b) \begin{align}&\partial {\phi}/\partial {r}=0,\quad r=b, \quad 0<z<\ell_{{n}}, \end{align}

(3.1c) \begin{align}& \partial {\phi}/\partial {z}=0,\quad z=0-,\quad b< r<a,\end{align}

(3.1d) \begin{align}&\partial {\phi}/\partial {r}=0,\quad r=a, \quad -\ell_{{w}}<z<0, \end{align}

(3.1e) \begin{align}& \phi(r,-\ell_{{w}})=0,\quad 0\leq r<a.\end{align}

This problem is similar to one considered in [1] and [10] (but note that we have interchanged a and b).

As in Section 2, the main quantity of interest is the flux,

\begin{equation*} {\mathcal J} =2\pi\int_0^{\rho(z)} \frac{\partial {\phi}}{\partial {z}}r\,{\rm d} r \quad\mbox{with}\quad \rho(z)=\left\{\begin{array}{lr} a, & -\ell_{{w}}<z<0,\\ \\[-7pt] b, & 0<z<\ell_{{n}}; \end{array}\right.\end{equation*}

${\mathcal J}$ does not depend on z. In particular, using (3.1c),

(3.2) \begin{equation} {\mathcal J}=2\pi\int_0^b \left.\frac{\partial {\phi}}{\partial {z}}\right|_{z=0} r\,{\rm d} r.\end{equation}

Applying Green’s theorem in the narrow tube to $\phi(r,z)$ and $z-\ell_{{n}}$ gives

\begin{equation*} 0=\pi b^2\phi_0-2\pi\int_0^b \phi(r,0+)r\,{\rm d} r -\ell_{{n}}{\mathcal J}.\end{equation*}

Similarly, applying Green’s theorem in the wide tube to $\phi(r,z)$ and $z+\ell_{{w}}$ gives

\begin{equation*} 0=2\pi\int_0^a \phi(r,0-)r\,{\rm d} r -\ell_{{w}}{\mathcal J}.\end{equation*}

Adding these two equations, noting that $\phi(r,0+)=\phi(r,0-)$ for $0\leq r<b$ , we obtain

(3.3) \begin{equation} (\ell_{{n}}+\ell_{{w}}){\mathcal J}=\pi b^2\phi_0 +2\pi\int_b^a \phi(r,0-)r\,{\rm d} r.\end{equation}

3.1 Boundary homogenisation

As in Section 2.1, suppose we approximate $\phi(r,z)$ by a function $\psi(z)$ , with $\psi''(z)=0$ . Imposing $\psi(\ell_{{n}})=\phi_0$ and $\psi(-\ell_{{w}})=0$ gives

\begin{equation*} \psi(z)=A(z+\ell_{{w}}),\quad -\ell_{{w}}<z<0 \quad\mbox{and}\quad \psi(z)=\phi_0 +B(z-\ell_{{n}}),\quad 0<z<\ell_{{n}},\end{equation*}

where A and B are arbitrary constants. The flux is $\pi a^2 A$ for $-\ell_{{w}}<z<0$ and $\pi b^2 B$ for $0<z<\ell_{{n}}$ . For these to match, $a^2A=b^2B$ . Thus, in terms of A,

(3.4) \begin{equation} \psi'(0-)=A, \quad \psi'(0+)=(a/b)^2A,\quad \psi(0-)=A\ell_{{w}},\quad \psi(0+)=\phi_0-(a/b)^2A\ell_{{n}}.\end{equation}

We need one more condition to determine A.

In [1, equation (2.3)], the following conditions are imposed at $z=0$ :

(3.5) \begin{equation} G_{{w}}^{\prime}(0-)=G_{{n}}^{\prime}(0+)=\kappa_{{n}} G_{{n}}(0+)-\kappa_{{w}} G_{{w}}(0-).\end{equation}

These two relations involve the functions $G_{{w}}$ and $G_{{n}}$ , and the parameters $\kappa_{{w}}$ and $\kappa_{{n}}$ .

The constants $\kappa_{{w}}$ and $\kappa_{{n}}$ are known as ‘boundary trapping rates’ [1]. Thus $\kappa_{{w}}$ is the trapping rate for diffusing particles approaching $z=0$ from the wide part of the structure (that is, from $z<0$ ) and $\kappa_{{n}}$ is the rate for particles arriving from the narrow part. These two rates are related, as we shall see.

The first of (3.5) is identified in [1] as flux conservation, so we put $G_{{w}}(z)=a^2\psi(z)$ for $-\ell_{{w}}<z<0$ and $G_{{n}}(z)=b^2\psi(z)$ for $0<z<\ell_{{n}}$ . The second of (3.5) then becomes

(3.6) \begin{equation} a^2\psi'(0-)=\kappa_{{n}} b^2\psi(0+)-\kappa_{{w}} a^2 \psi(0-).\end{equation}

In [1, equation (2.5)], it is argued that the boundary trapping rates satisfy

(3.7) \begin{equation} \kappa_{{n}} = (a/b)^2\kappa_{{w}} =\kappa,\end{equation}

say. (The authors appeal to ‘the condition of detailed balance’ but an elementary argument suffices. If we replace the boundary condition (3.1e) by $\phi(r,-\ell_{{w}})=\phi_0$ , the exact solution is $\phi(r,z)=\phi_0$ with zero flux. The corresponding one-dimensional solution is $\psi(z)=\phi_0$ , but this solution will only satisfy (3.6) if (3.7) holds.) Using (3.7) reduces (3.6) to

(3.8) \begin{equation} a^2\psi'(0-)=b^2\kappa \left\{ \psi(0+)-\psi(0-)\right\}.\end{equation}

We notice that the right-hand side of (3.8) contains the discontinuity in $\psi$ across the junction at $z=0$ whereas $\phi(r,0+)=\phi(r,0-)$ for $0\leq r<b$ .

Substituting from (3.4) in (3.8), we find that $A=\kappa b^2\phi_0/(a^2+\kappa\Omega)$ and

(3.9) \begin{equation} {\mathcal J}=\pi a^2 A=\frac{\pi a^2b^2\phi_0}{a^2\kappa^{-1}+\Omega} \quad\mbox{with}\quad \Omega=a^2\ell_{{n}}+b^2\ell_{{w}}.\end{equation}

It is known [1] that $\kappa\to\infty$ as $b\to a$ , and so the correct result is obtained in this limit.

3.2 Two semi-infinite tubes

Returning to problem (3.1), suppose we let $\ell_{{w}}\to\infty$ and $\ell_{{n}}\to\infty$ , and consider uniform flow along two coaxial semi-infinite tubes, joined together with a rigid annulus at $z=0$ . For this problem, the boundary conditions (3.1a) and (3.1e) are replaced by far-field conditions,

\begin{equation*} \phi(r,z)=\left\{\begin{array}{l@{\quad}l@{\quad}l} U_{{n}}(z+C)+o(1)\quad &\mbox{as $z\to\infty$,} &0\leq r<b,\\ \\[-6pt] Uz+o(1)\quad& \mbox{as $z\to-\infty$,} & 0\leq r<a. \end{array}\right.\end{equation*}

Flux conservation shows that the constants U and $U_{{n}}$ are related by $\pi a^2U=\pi b^2 U_{{n}}$ . The constant C is unknown; it is the blockage coefficient. As all boundary conditions are Neumann conditions, we are free to exclude an additive constant in the behaviour of $\phi(r,z)$ as $z\to -\infty$ .

If we apply Green’s theorem to $\phi(r,z)$ and z in the narrow tube, we obtain

\begin{equation*} U_{{n}} C\pi b^2=2\pi\int_0^b \phi(r,0+)r\,{\rm d} r.\end{equation*}

If we do the same in the wide tube, we obtain

\begin{equation*} 2\pi\int_0^a \phi(r,0-)r\,{\rm d} r =0.\end{equation*}

As $U_{{n}} b^2=Ua^2$ and $\phi(r,0+)=\phi(r,0-)$ for $0\leq r<b$ , we obtain

(3.10) \begin{equation} C= -\frac{2}{Ua^2}\int_b^a \phi(r,0-)r\,{\rm d} r.\end{equation}

Our goal is to obtain an approximation to C; (3.10) shows that we are mainly interested in $\phi$ in the wide tube. (Of course, a full solution would require matching to $\phi$ in the narrow tube through the aperture at $z=0$ ; see [20] for some details.)

As the flow is axisymmetric, write $\phi$ in the wide tube as

(3.11) \begin{equation} \phi(r,z)=Uz +Ua\sum_{n=1}^\infty c_n J_0(\lambda_nr)\,{\rm e}^{\lambda_nz}, \quad z<0,\quad 0\leq r<a,\end{equation}

where $J_1(\lambda_na)=0$ , as before. Differentiating,

\begin{equation*} U+U\sum_{n=1}^\infty c_n(\lambda_na) J_0(\lambda_nr)=\left\{ \begin{array}{c@{\quad}l} V(r), & 0\leq r<b\\ \\[-7pt] 0, & b<r<a, \end{array}\right.\end{equation*}

where V(r) is the unknown speed in the aperture, $V(r) =\partial {\phi}/\partial {z}$ at $z=0$ , $0\leq r<b$ and we have used (3.12c). Orthogonality gives

(3.12) \begin{equation} \pi a^2U=2\pi\int_0^bV(r)r\,{\rm d} r\end{equation}

(which is flux conservation again) and

\begin{equation*} c_n =\frac{2}{U\lambda_na^3\,J_0^2(\lambda_na)}\int_0^b V(r)J_0(\lambda_nr)r\,{\rm d} r.\end{equation*}

Also, substituting (3.11) in (3.10) gives

(3.13) \begin{align}C&= -\frac{2}{a}\sum_{n=1}^\infty c_n \int_b^aJ_0(\lambda_nr)r\,{\rm d} r= \frac{2b}{a}\sum_{n=1}^\infty \frac{c_n}{\lambda_n} J_1(\lambda_nb) \nonumber\\[3pt] &=\frac{4b}{Ua^2}\sum_{n=1}^\infty \frac{J_1(\lambda_nb)}{(\lambda_na)^2 J_0^2(\lambda_na)}\int_0^b V(r)J_0(\lambda_nr)r\,{\rm d} r, \end{align}

using (2.12). The formula (3.13) is exact, but we do not know V(r).

To proceed, we approximate and try $V(r)=U_q(b^2-r^2)^{-q}$ , with $0\leq q<1$ ; conservation gives $U_q=(1-q)a^2b^{2q-2}U$ . A local analysis near the corner at $(r,z)=(b,0)$ gives $q=\frac{1}{3}$ [12, Section 63]. Note that $q=0$ would mean that V(r) has been approximated by a constant whereas $q=\frac{1}{2}$ would be appropriate when two identical semi-infinite tubes are joined with a rigid iris (a thin screen with a hole of radius b). Then (3.13) gives

(3.14) \begin{align}C \simeq C_{{\rm tube}}^{{\rm app}}&= \frac{4bU_q}{Ua^2}\sum_{n=1}^\infty \frac{J_1(\lambda_nb)}{(\lambda_na)^2 J_0^2(\lambda_na)}\int_0^b (b^2-r^2)^{-q} J_0(\lambda_nr)r\,{\rm d} r \nonumber\\&=2^{2-q} b \,\Gamma(2-q)\,(b/a)^{q-1}\,{\mathbb S}_q(b/a),\end{align}

where

(3.15) \begin{equation} {\mathbb S}_q(X)= \sum_{n=1}^\infty \frac{J_1(\lambda_naX)\,J_{1-q}(\lambda_naX)}{(\lambda_na)^{3-q}\,J_0^2(\lambda_na)} =\sum_{n=1}^\infty \frac{J_1(j_nX)\,J_{1-q}(j_nX)}{j_n^{3-q}J_0^2(j_n)},\end{equation}

$j_n=\lambda_na$ and we have used [7, 6.567.1]

(3.16) \begin{equation} \int_0^b (b^2-r^2)^{-q} J_0(\lambda_nr)r\,{\rm d} r =\frac{b^{2-2q}\,\Gamma(1-q)}{2^q (\lambda_nb)^{1-q}}\,J_{1-q}(\lambda_nb).\end{equation}

Note that the series ${\mathbb S}_0$ appeared earlier; see (2.17). Asymptotic approximations to C for $b/a\ll 1$ are developed in Appendix A.

3.3 Junction problem: Estimating $\kappa$

Let us return to the junction problem (3.1).

3.3.1 The wide tube

In the wide tube, we can write

(3.17) \begin{equation} \frac{\phi(r,z)}{\phi_0}=c_0\left(\frac{z}{\ell_{{w}}}+1\right) +\sum_{m=1}^\infty c_m J_0(\lambda_mr)\sinh{\{\lambda_m(z+\ell_{{w}})\}}\end{equation}

for $-\ell_{{w}}<z<0$ and $0\leq r<a$ , where $J_1(\lambda_ma)=0$ . This representation satisfies (3.1d) and (3.1e). At $z=0$ , we have

(3.18) \begin{equation} c_0+\sum_{m=1}^\infty c_m J_0(\lambda_mr)\sinh{\lambda_m\ell_{{w}}}=\frac{\varphi(r)}{\phi_0}, \quad 0\leq r<b,\end{equation}

where $\varphi(r)=\phi(r,0)$ , $0\leq r<b$ ; $\varphi(r)$ is the unknown potential in the aperture. Differentiating (3.17),

(3.19) \begin{equation} c_0+\sum_{m=1}^\infty c_m J_0(\lambda_mr)(\lambda_m\ell_{{w}})\cosh{\lambda_m\ell_{{w}}} =\left\{\begin{array}{cl} (\ell_{{w}}/\phi_0) V(r), &0\leq r<b, \\ \\[-7pt] 0, & b<r<a, \end{array}\right.\end{equation}

where $V(r)=\partial {\phi}/\partial {z}$ at $z=0$ , $0\leq r<b$ and we have used (3.1c). Orthogonality gives

(3.20) \begin{equation} c_0=\frac{2\ell_{{w}}}{\phi_0a^2}\int_0^b V(r)r\,{\rm d} r =\frac{\ell_{{w}}}{\phi_0\pi a^2}\,{\mathcal J}\end{equation}

using (3.2) and, using (2.15),

\begin{equation*}c_m \lambda_m\ell_{{w}}\cosh{\lambda_m\ell_{{w}}}=\frac{2\ell_{{w}}}{\phi_0a^2 J_0^2(\lambda_ma) } \int_0^b V(r) J_0(\lambda_mr) r\,{\rm d} r.\end{equation*}

If we eliminate $c_m$ from (3.17) and assume that V(r) is constant, we obtain [11, equation (3)]. Instead of doing that, we eliminate $c_m$ from (3.18), giving an equation relating two unknown functions, V and $\varphi$ ,

(3.21) \begin{equation} \frac{2\ell_{{w}}}{a^2}\int_0^b V(s) \bigg\{ 1 + \frac{a}{\ell_{{w}}}\sum_{m=1}^\infty \frac{J_0(\lambda_mr) J_0(\lambda_ms)}{ (\lambda_ma) J_0^2(\lambda_ma)}\tanh{\lambda_m\ell_{{w}}} \bigg\} s\,{\rm d} s =\varphi(r), \quad 0\leq r<b.\end{equation}

3.3.2 The narrow tube

Let us make a similar calculation for the narrow tube, starting with

(3.22) \begin{equation} \frac{\phi(r,z)}{\phi_0}=1+d_0\left(\frac{z}{\ell_{{n}}}-1\right) +\sum_{m=1}^\infty d_m J_0(\mu_m r)\sinh{\{\mu_m(z-\ell_{{n}})\}},\end{equation}

for $0<z<\ell_{{n}}$ and $0\leq r<b$ , where $\mu_m$ are positive solutions of $J_1(\mu_mb)=0$ and the coefficients $d_m$ are to be determined. This representation for $\phi$ ensures that (3.1a) and (3.1b) are satisfied. At $z=0$ ,

(3.23) \begin{equation} 1-d_0-\sum_{m=1}^\infty d_m J_0(\mu_mr)\sinh{\mu_m \ell_{{n}}}=\frac{\varphi(r)}{\phi_0}, \quad 0\leq r<b.\end{equation}

Differentiating (3.22),

(3.24) \begin{equation} d_0 +\sum_{m=1}^\infty d_m J_0(\mu_mr) (\mu_m\ell_{{n}})\cosh{\mu_m\ell_{{n}}} =(\ell_{{n}}/\phi_0) V(r), \quad 0\leq r<b.\end{equation}

From this equation, we obtain

(3.25) \begin{equation} d_0=\frac{2 \ell_{{n}}}{\phi_0b^2} \int_0^b V(r)r\,{\rm d} r =\frac{\ell_{{n}}}{\phi_0\pi b^2}\,{\mathcal J}\end{equation}

and, using (2.15),

\begin{equation*}d_m \mu_m\ell_{{n}}\cosh{\mu_m\ell_{{n}}}=\frac{2\ell_{{n}}}{\phi_0b^2 J_0^2(\mu_mb)}\int_0^b V(r) J_0(\mu_mr)r\,{\rm d} r.\end{equation*}

Eliminating $d_m$ from (3.23),

(3.26) \begin{equation} -\frac{2 \ell_{{n}}}{b^2} \int_0^b V(s) \bigg\{1 + \frac{b}{\ell_{{n}}} \sum_{m=1}^\infty\frac{ J_0(\mu_mr) J_0(\mu_ms) }{(\mu_mb) J_0^2(\mu_mb)} \tanh{\mu_m\ell_{{n}}} \bigg\}s\,{\rm d} s =\varphi(r)-\phi_0 \quad 0\leq r<b.\end{equation}

Chapman and Parker [3] eliminate $\varphi$ between (3.18) and (3.23) and V between (3.19) and (3.24), giving a set of dual series equations for $c_m$ and $d_m$ . They solve these equations numerically, with a focus on computing $c_0$ , which is related to the flux by (3.20). See also [20, Section III.B].

3.3.3 An integral equation for V

Eliminating $\varphi$ between (3.21) and (3.26), we obtain

(3.27) \begin{align}& \frac{2\ell_{{w}}}{a^2}\int_0^b V(s) \bigg\{ 1 + \frac{a}{\ell_{{w}}}\sum_{m=1}^\infty \frac{J_0(\lambda_mr) J_0(\lambda_ms)}{ (\lambda_ma) J_0^2(\lambda_ma)}\tanh{\lambda_m\ell_{{w}}} \bigg\} s\,{\rm d} s \nonumber\\[3pt] & \quad\mbox{} +\frac{2 \ell_{{n}}}{b^2} \int_0^b V(s) \bigg\{1 + \frac{b}{\ell_{{n}}} \sum_{m=1}^\infty\frac{ J_0(\mu_mr) J_0(\mu_ms) }{(\mu_mb) J_0^2(\mu_mb)} \tanh{\mu_m\ell_{{n}}} \bigg\}s\,{\rm d} s =\phi_0,\end{align}

which is an exact integral equation for V(r), $0\leq r<b$ . It is unlikely that analytical solutions can be found, so we seek approximations.

3.3.4 Approximation for long tubes

When $\ell_{{w}}/a$ and $\ell_{{n}}/b$ are large, we can replace both hyperbolic tangents in (3.27) by 1. Moreover, at leading order, we can discard the infinite series because they are multiplied by $a/\ell_{{w}}$ or $b/\ell_{{n}}$ . Doing this leaves the simple approximation

(3.28) \begin{equation}\left(\frac{2\ell_{{w}}}{a^2}+\frac{2\ell_{{n}}}{b^2}\right)\int_0^bV(s)s\,{\rm d} s=\phi_0\end{equation}

which gives

(3.29) \begin{equation} {\mathcal J}\simeq \pi a^2b^2\phi_0/\Omega \quad\mbox{where}\quad \Omega=a^2\ell_{{n}}+b^2\ell_{{w}}\end{equation}

appeared in (3.9). Also, comparing (3.29) with (3.20) gives

(3.30) \begin{equation} c_0\simeq b^2\ell_{{w}}/\Omega.\end{equation}

To improve our approximation, we proceed as in Section 2.3.2 and use an iterative approach. Suppose that $V=V_0\ell_{{w}}^{-1}+V_1\ell_{{w}}^{-2}+\cdots$ , where $V_0$ is chosen so that $V(s)=V_0(s)\ell_{{w}}^{-1}$ satisfies (3.28). If we choose $V_0(r)=W_q(b^2-r^2)^{-q}$ , we obtain $W_q=\phi_0\ell_{{w}} a^2b^{2q}(1-q)/\Omega$ . Hence, we find that $V_1$ satisfies

\begin{align*}& \frac{2\ell_{{w}}}{a^2}\int_0^b \bigg\{ \frac{V_1(s)}{\ell_{{w}}^2} + \frac{aV_0(s)}{\ell_{{w}}^2}\sum_{m=1}^\infty \frac{J_0(\lambda_mr) J_0(\lambda_ms)}{ (\lambda_ma) J_0^2(\lambda_ma)} \bigg\} s\,{\rm d} s \\[3pt] &\quad\mbox{}+ \frac{2 \ell_{{n}}}{b^2} \int_0^b \bigg\{\frac{V_1(s)}{\ell_{{w}}^2} + \frac{bV_0(s)}{\ell_{{n}}\ell_{{w}}} \sum_{m=1}^\infty\frac{ J_0(\mu_mr) J_0(\mu_ms) }{(\mu_mb) J_0^2(\mu_mb)} \bigg\}s\,{\rm d} s =0,\quad 0\leq r<b.\end{align*}

Using (3.16), we have

\begin{equation*} \int_0^b V_0(s)\,J_0( \lambda s)s\,{\rm d} s=\frac{\phi_0\ell_{{w}} a^2b^2 \,\Gamma(2-q)}{2^q\,\Omega\,(\lambda b)^{1-q}}\,J_{1-q}(\lambda b)\end{equation*}

whence

\begin{align*} \frac{2}{\ell_{{w}}^2}\left(\frac{\ell_{{w}}}{a^2}+\frac{\ell_{{n}}}{b^2}\right)\int_0^b V_1(s)s\,{\rm d} s&=-\frac{\phi_0a^2b (b/a)^q\Gamma(2-q)}{2^{q-1}\,\Omega} \sum_{m=1}^\infty \frac{J_0(\lambda_mr)\,J_{1-q}(\lambda_mb)}{(\lambda_ma)^{2-q}J_0^2(\lambda_ma)}\\[3pt] &\quad\mbox{} - \frac{\phi_0a^2b \,\Gamma(2-q)}{2^{q-1}\,\Omega} \sum_{m=1}^\infty \frac{J_0(\mu_mr)\,J_{1-q}(\mu_mb)}{(\mu_ma)^{2-q}J_0^2(\mu_ma)}.\end{align*}

But, as the right-hand side depends on r, we average over the disc: multiply by $2r/b^2$ and integrate between $r=0$ and $r=b$ , using

\begin{equation*} \int_0^b J_0(\lambda_mr)r\,{\rm d} r=\frac{b}{\lambda_m}J_1(\lambda_mb) \quad\mbox{and}\quad \int_0^b J_0(\mu_mr)r\,{\rm d} r =0.\end{equation*}

Doing this gives

\begin{equation*} \frac{2}{\ell_{{w}}^2}\left(\frac{\ell_{{w}}}{a^2}+\frac{\ell_{{n}}}{b^2}\right)\int_0^b V_1(s)s\,{\rm d} s =-\frac{\phi_0a^3 (b/a)^q\Gamma(2-q)}{2^{q-2}\,\Omega} \,{\mathbb S}_q(b/a),\end{equation*}

where ${\mathbb S}_q(X)$ is defined by (3.15). For the flux, we obtain

\begin{align*} {\mathcal J}&=2\pi\int_0^bV(s)s\,{\rm d} s=2\pi\int_0^b \left(\frac{V_0}{\ell_{{w}}}+\frac{V_1}{\ell_{{w}}^2}\right)s\,{\rm d} s\\[3pt] & =\frac{\pi a^2 b^2\phi_0}{\Omega}\left(1 -\frac{a^3(b/a)^q}{2^{q-2}\Omega}\,\Gamma(2-q)\,{\mathbb S}_q(b/a) \right),\end{align*}

where $\Omega$ is defined by (3.29). The sum appearing here also appears in a comparable approximation for the blockage coefficient C for uniform flow along two coaxial tubes joined together, $C_{{\rm tube}}^{{\rm app}}$ , defined by (3.14). Using this, we obtain

(3.31) \begin{equation} {\mathcal J} =\frac{\pi a^2\phi_0}{\ell_{{w}}}\,c_0\simeq \frac{\pi a^2 b^2\phi_0}{\Omega}\left(1 -\frac{a^2C_{{\rm tube}}^{{\rm app}}}{\Omega} \right).\end{equation}

3.3.5 Estimating $\kappa$

We have two estimates for ${\mathcal J}$ , namely (3.9) and (3.31). Equating these gives

\begin{equation*}\frac{1}{a^2\kappa^{-1}+\Omega} =\frac{1}{\Omega}\left(1 -\frac{a^2C_{{\rm tube}}^{{\rm app}}}{\Omega} \right).\end{equation*}

Solving for $\kappa$ , we obtain

(3.32) \begin{equation} \kappa=\frac{1}{C_{{\rm tube}}^{{\rm app}}}-\frac{a^2}{a^2\ell_{{n}}+b^2\ell_{{w}}}.\end{equation}

Again, we see the slow decay of $\kappa$ with increasing tube lengths.

If we combine the leading-order estimate for $\kappa$ with (3.7), we obtain

\begin{equation*} \kappa_{{w}} = \frac{b^2}{a^2}\kappa \simeq \frac{b^2}{a^2C_{{\rm tube}}^{{\rm app}}}. \end{equation*}

For small $b/a$ , we can use the asymptotic approximation (A3), giving

(3.33) \begin{equation} \kappa_{{w}}\simeq \frac{bA}{a^2}\quad\mbox{with}\quad A=\left(\frac{3}{2}-q\right)\!\left(\frac{\Gamma(\frac{3}{2}-q)}{\Gamma(2-q)}\right)^2,\end{equation}

where it was assumed that V(r) is proportional to $(b^2-r^2)^{-q}$ with $0\leq q<1$ . The proper choice for q is $q=\frac{1}{3}$ ; this gives $A\simeq 1.23$ .

If we take $q=0$ (so that V(r) is approximated by a constant, the average flux into the narrow pipe), we obtain $A=\frac{3}{8}\pi\simeq 1.18$ .

If we take $q=\frac{1}{2}$ , we obtain $A=4/\pi\simeq 1.27$ . This is the limiting value in [1, equation (2.4a)]. The inverse square-root behaviour at $r=b$ implies that the underlying flow problem has been replaced by the problem of flow through a hole in a thin rigid screen, leading to an error that can be estimated using (3.33). (For details and references on this ‘iris problem’, see [14, 15].) Qualitatively, this is a mistake because it ignores the presence of the narrow tube; however, quantitatively, the values of A for $q=\frac{1}{3}$ ( $A\simeq 1.23$ ) and $q=\frac{1}{2}$ ( $A\simeq 1.27$ ) are close.

Of course, the assumption that V(r) is proportional to $(b^2-r^2)^{-q}$ is itself an approximation: in principle, the exact solution for V could be found by solving the integral equation (3.27), but this seems to be difficult analytically without further approximations. Numerical solutions could be sought, but doing this would be expensive (and it would be outside the scope of this paper).

3.4 Junction problem: The one-dimensional approximation $\boldsymbol\psi$

Let us construct $\psi$ from $\phi$ , using (1.1), which reduces to

(3.34a) \begin{align} \psi(z)&=\frac{2}{a^2}\int_0^a \phi(r,z)r\,{\rm d} r=\phi_0c_0\left(\frac{z}{\ell_{{w}}}+1\right), \qquad -\ell_{{w}}<z<0, \end{align}

(3.34b) \begin{align} \psi(z)&=\frac{2}{b^2}\int_0^b \phi(r,z)r\,{\rm d} r=\phi_0\left\{1+d_0\left(\frac{z}{\ell_{{n}}}-1\right) \right\}, \quad 0<z<\ell_{{n}},\end{align}

after use of (3.17) and (3.22). These satisfy $\psi(-\ell_{{w}})=0$ and $\psi(\ell_{{n}})=\phi_0$ . Flux continuity, $a^2\psi'(0-)=b^2\psi'(0+)$ gives $a^2c_0/\ell_{{w}}=b^2d_0/\ell_{{n}}$ , which is seen to hold; compare (3.20) and (3.25). Then, applying the condition (3.8) at $z=0$ gives

\begin{equation*} a^2c_0=\kappa b^2\ell_{{w}}(1-c_0-d_0)=\kappa(b^2\ell_{{w}}-\Omega c_0),\end{equation*}

which gives an exact formula for $\kappa$ ,

(3.35) \begin{equation} \kappa=\frac{a^2c_0}{b^2\ell_{{w}}-\Omega c_0}.\end{equation}

Now, the leading-order estimate for $c_0$ , (3.30), is $c_0\simeq b^2\ell_{{w}}/\Omega$ , and so (3.35) does not give a useful estimate of $\kappa$ . However, if we use the refined estimate for $c_0$ , (3.31), we recover (3.32).

As ${\mathcal J}=\pi a^2\phi_0c_0/\ell_{{w}}=\pi b^2\phi_0d_0/\ell_{{n}}$ , we can write (3.34) as

\begin{align*}& \pi a^2\psi(z)={\mathcal J}(z+\ell_{{w}}),\qquad -\ell_{{w}}<z<0,\\[3pt] & \pi b^2\psi(z)={\mathcal J}(z-\ell_{{n}} +\ell_{{n}}/d_0),\quad 0<z<\ell_{{n}}.\end{align*}

Kalinay and Percus [10, equation (39)] obtained similar formulas. However, their analysis of (3) does not take account of finite tube lengths and does not gives estimates of $\kappa$ .

4 Discussion

We have given a detailed study or two problems involving what are nominally one-dimensional geometries: long tubes with mixed boundary conditions at one end (the ‘finite-tube problem’ discussed in Section 2) or with an abrupt change in cross-section (the ‘junction problem’ discussed in Section 3). Such problems are often tackled using one-dimensional models coupled with certain effective boundary conditions; generically, these conditions contain a parameter $\kappa$ . The same problems can also be tackled using eigenfunction expansions (separation of variables). This approach is more complicated but, in principle, it is exact: it can accommodate all lateral variations of the solution. Comparing the two approaches has given us some insight into $\kappa$ and its properties. For long tubes, we might expect to see connections with the blockage coefficient C; this quantity is uniquely defined by solving a related potential flow problem. Indeed, we find that, generically, $\kappa$ is proportional to $1/C$ . But we also found that $\kappa$ depends on the length of the tube $\ell$ ; it approaches its limiting value slowly (as $\ell^{-1}$ ) as $\ell$ increases. We also saw that moving the location of the effective boundary condition (from $z=z_0$ to $z=z_1$ ) also affects $\kappa$ ; see the text around (2.27). These properties mean that we cannot view $\kappa$ as being an intrinsic property of the boundaries or interfaces being modelled.

A further difficulty comes when $\kappa$ is viewed as a boundary property. For example, $\kappa$ has been related to ‘the Hill formula [8] for the flux to an isolated disk on a reflecting wall’ [2, p. 2]. In detail, Hill [8, p. 4919] solved $\nabla^2\phi=0$ in the half-space $z>0$ with a rigid/reflective boundary at $z=0$ apart from a circular ‘window’ on which he imposed the boundary condition $\phi=0$ . This implies that $\phi$ can be extended into the region $z<0$ as an odd function of z, $\phi(r,-z)=-\phi(r,z)$ . In other words, the window is regarded as a circular aperture in a thin rigid plane (leading to a problem that Hill could solve using oblate spheroidal coordinates). But the window could be the mouth of a tube extending into $z<0$ , or maybe a cone; in situations such as these, we cannot impose $\phi=0$ on the window, we have to impose continuity conditions between the two regions, just as we did for the finite-tube problem (Section 2) and for the junction problem (Section 3). Thus, the geometry on both sides of the reflecting wall should be taken into account: the conceptual replacement of the entry cross-section of the narrow tube by an absorbing disc can lead to errors. Similar errors are to be expected when one-dimensional models are used for more complicated geometrical configurations. Quantifying these errors is difficult, because they depend on the parameters and geometry of the specific problem of interest: care is warranted.