## Abstract

We present a linear stability analysis of viscous, pressure-driven flow in an annular pipe with a moving core. The core moves with constant speed along its axis and rotates with constant angular frequency. At both boundaries the fluid obeys slip boundary conditions.

Our numerical results highlight the important role of asymmetric perturbations in the instability of the system. Boundary slip is seen to be stabilizing, as has been observed for other types of parallel flows. We show that motion of the solid core with respect to the outer cylinder causes the neutral curve to split into two distinct curves. In particular, we show that rotation has a surprising destabilizing effect.

## 1. Introduction

The papers of Lüscher *et al*. (1995) and Frei *et al*. (2000) have succeeded in drawing attention back to an important problem from the field of hydrodynamic stability. Concerned with how to introduce a porous medium into a patient in a minimally invasive way, Frei *et al*. (2000) provide a detailed description of the modern technique of *thread injection*, in which a hypodermic needle is connected to a syringe of fluid, which also contains a spool of surgical thread. One end of the thread is introduced into the needle, which in turn is inserted into the patient. A pressure gradient is applied by the syringe, and the thread is drawn through the needle with the fluid into the patient, the speed of the thread being controlled by a motor connected to the spool. The *thread-annular problem* is the investigation into the stability of this type of fluid system and is the motivation for this work.

Research into the thread-annular problem can be seen to have started with the work of Mott & Joseph (1968), who presented a linear stability analysis for a viscous fluid flowing along a constant pressure gradient through a pipe with a solid core along its axis. Letting the pipe have inner radius *R*, the core have radius *δ* and setting *η*=*δ*/*R*, Mott & Joseph investigated the instability of the system to axisymmetric perturbations, tracking the critical Reynolds number guaranteeing instability as a function of *η*. Their results showed the critical Reynolds number varying monotonically with *η* from a finite value at *η*=1 (corresponding to Poiseuille flow in the plane), approaching infinity as *η*→0 (corresponding to Hagen–Poiseuille pipe flow, see the book by Drazin & Reid (1981) for a discussion). This was verified later in the work of Strumolo (1982).

More recent have been the contributions of Walton (2003, 2004, 2005), whose model allows the inner cylinder to move axially with constant velocity *W* in the direction of the pressure gradient, and for algebraic neatness uses a non-dimensionalization different from that of Mott & Joseph. Walton showed that multiple neutral curves can coexist when *W*≠0, making the problem of tracking the critical Reynolds number much more difficult.

We bring this model closer to the problem considered by Frei *et al*. (2000) by including two further types of motion. The first is we allow the solid core to rotate about the axis with constant angular frequency *f* and show that this introduces another component into the steady-state solution, and that increasing the rate of rotation has a surprising, destabilizing effect on the stability of the system.

Finally, we introduce slip boundary conditions, which we allow, since we consider our model to be valid on the interior of a hypodermic needle. There is increasing evidence to suggest that in some situations the no-slip boundary conditions for the Navier–Stokes equations may be inappropriate. At a fluid–solid interface, there is some motion of the fluid relative to the solid surface, which has been observed experimentally (see Craig *et al*. 2001, 2003; Zhu & Granick 2002; Tabeling 2004; Tabeling & Joseph 2005). Therefore, there has recently been interest in incorporating slip boundary effects into models of fluid flow, for instance Jie *et al*. (2000), Yu & Ameel (2001), De Socio & Marino (2002) and Raju & Roy (2004).

Central to the slip boundary condition is the concept of the ‘slip length’ *λ*, which is of the order of the distance beyond the boundary at which the fluid velocity extrapolates to zero. Therefore, boundary slip is significant in geometries where the characteristic length is small, approaching the size of *λ*. Hypodermic needles in common use can have inner diameters as small as 76.2×10^{−6} m, and so we believe boundary slip to be an important feature of the thread-annular model.

We begin by describing our model mathematically, deriving non-traditional boundary conditions and obtaining a steady-state solution, before performing a linear stability analysis. We then present the results of our detailed, computational investigation of the effects of each parameter on the stability of the system.

## 2. Governing equations

Adopting the cylindrical polar coordinate system (*r*, *θ*, *z*), we consider a viscous, incompressible fluid on the annular domain *Ω*_{R}={(*r*, *θ*, *z*)|*r*∈[*δ*,*R*],*θ*∈[0,2*π*),*z*∈(−∞, ∞)}. At time *t*∈[−0,∞) a fluid particle has velocity , which satisfies the Navier–Stokes equations on *Ω*_{R},(2.1)(2.2)(2.3)(2.4)where *p*(*r*, *θ*, *z*; *t*) is the pressure field, and *ρ* and *ν* are the constant density and kinematic viscosity, respectively.

The thread-annular model is introduced through imposing appropriate conditions on the fluid at the surfaces *r*=*δ* and *r*=*R*. We suppose that the thread moves with constant velocity *W*≥0 along the *z*-axis, and rotates with angular frequency *f*≥0, such that points on the surface of the thread move with velocity . The surface *r*=*R* is motionless.

The slip boundary condition (see Spille *et al*. 2000; Lauga & Cossu 2005; Webber 2006; Webber & Straughan 2006) is to impose a zero-flux condition at the boundaries, i.e. *u*|_{r=δ}=*u*|_{r=R}=0, and set the tangential components *v* and *w* to be proportional to the corresponding components of the shear stress ±*ϵ*_{ir}, i.e.(2.5)(2.6)where *λ*_{R}≥0 is the constant of proportionality, or slip length, at *r*=*R* and *λ*_{δ} is defined similarly. Once we also allow for the motion of the thread, we arrive at the boundary conditions for our model or the thread-annular problem,(2.7)

(2.8)

### (a) Steady-state solution

We consider a time-independent solution of equations (2.1)–(2.8) of the form(2.9)arising from the shearing effect of the thread's motion and a constant pressure gradient . Substituting these forms into equation (2.4), we obtain(2.10)where the prime denotes differentiation with respect to *r*. Since we obtain . We are left to solve the system(2.11)(2.12)(2.13)By differentiating equation (2.11) with respect to *θ*, we find that is constant. Therefore we impose , else would vary with *θ* as a sawtooth wave, discontinuous at *θ*=2*π*.

Finally, reconciling the boundary conditions at *r*=*δ* and *R*, we obtain the dimensional steady-state solution(2.14)

(2.15)

(2.16)

### (b) Non-dimensionalization

For algebraic neatness, we non-dimensionalize using the scalings of Walton (2003, 2004, 2005) rather than employing the scalings common to the treatment of parallel flow problems. Denoting non-dimensional quantities by the notation ‘^{*}’, the scalings we use arewhere *Re* is our Reynolds number.

Following non-dimensionalization and upon dropping the ‘^{*}’ notation, our equations of motion (2.1)–(2.4) become the familiar non-dimensional Navier–Stokes equations(2.17)(2.18)(2.19)(2.20)valid on the domain *Ω*={(*r*, *θ*, *z*)|*r*∈[*δ*,1],*θ*∈[0,2*π*),*z*∈(−∞, ∞)}, which yield the steady-state solution(2.21)

(2.22)

(2.23)

### (c) Instability equations

Our aim is to investigate the instability of the steady-state to three-dimensional disturbances . We begin by perturbing the steady-state , , and substituting this perturbed form into equations (2.17)–(2.20). Upon removing terms nonlinear in ** u**, we obtain the linearized perturbation equations as follows:(2.24)(2.25)(2.26)(2.27)where the components of

**satisfy the boundary conditions(2.28)**

*u*(2.29)

Rather than performing the stability analysis in the primitive variables (*u*, *v*, *w*, *p*), we employ an Orr–Sommerfeld technique. Although this sacrifices simplicity of coding and incurs the presence of spurious eigenvalues, the Orr–Sommerfeld approach has been seen to make convergence faster (see Kerswell & Davey 1996) and spurious eigenvalues may easily be identified and filtered out of our algorithm.

The Orr–Sommerfeld system we arrive at, and the way in which we obtain it, depends on whether we are interested in symmetric or asymmetric perturbations. This is due to the way *n* appears from *∂*/*∂**θ* terms once we make our modal substitution (see below), therefore we must derive both symmetric and asymmetric Orr–Sommerfeld equations. We outline the treatment for the asymmetric case, but the case with symmetric perturbations follows similarly. See the book by Drazin & Reid (1981) for more examples.

To remove the pressure terms, we perform the operationsEquation (2.27) is then used to eliminate the *w* component, leaving two equations in *u* and *v*,(2.30)

(2.31)where, by transferring the boundary conditions upon *w* onto *u*, *u* and *v* can be shown to satisfy(2.32)

(2.33)

We suppose that *u*(*r*, *θ*, *z*;*t*) can be expressed as an expansion of Fourier modes, specifically(2.34)with a similar expression for *v*(*r*, *θ*, *z*;*t*), for positive wavenumbers *a*_{m} and complex growth factors *c*_{m}. However, since unbounded growth of a single mode is sufficient for instability, we need to only consider travelling wave solutions of the form(2.35)

Substitution of these forms into equations (2.30) and (2.31) leaves us with a system of two coupled differential equations in *r* with eigenvalue *c* and eigenfunction {*u*(*r*),*v*(*r*)}.

## 3. Numerical solution

We solve system (2.30) and (2.31) using the Chebyshev-tau method of Orszag (1971), which we shall briefly describe below. The method exploits the orthogonality of Chebyshev polynomials with respect to the weighted scalar product on *L*^{2}(−1, 1),(3.1)where *T*_{k}(*x*) is the Chebyshev polynomial of order *k*≥0.

We begin, therefore, by mapping the problem from *r*∈[*δ*,1] to *x*∈[−1,1] and expanding the functions *u*(*x*), *v*(*x*) and their derivatives as infinite series of Chebyshev polynomials, e.g. , , and so on. We then substitute these forms into (2.30) and (2.31); note that the coefficients , can be expressed in terms of the coefficients *u*_{m} (see Orszag), and the same is the case with *v*_{m}, before truncating these expansions at the *k*+1th term.

This gives us a problem in 2*k*+2 unknowns. By multiplying our truncated forms in (2.30) and (2.31) with the polynomials *T*_{m}(*x*) for *m*=0, …, *k*−3, and in each case taking the scalar product (3.1), we obtain 2*k*−4 linear equations. The system is closed using the six boundary conditions in (2.32) and (2.33). Ultimately, we arrive at the generalized eigenvalue problem(3.2)where *x*^{T}=(*u*_{0}, *u*_{1}, …, *u*_{k}, *v*_{0}, *v*_{1}, …, *v*_{k}); *A* and *B* are square matrices with complex entries; and *c* is the complex growth factor from (2.35), which appears as a result of the *∂*/*∂**t* terms in (2.30) and (2.31). We then solve this using the QZ algorithm of Moler & Stewart (1973) via the Numerical Algorithms Group (NAG) routine F02GJF.

We have only outlined the basics of the Chebyshev-tau method. As well as the original paper of Orszag (1971), the reader is referred to the work of Dongrarra *et al*. (1996) and the book by Straughan (2004) for practical details of implementing the method. Under the nomenclature of Dongarra *et al*. we use the *D*^{2} implementation of the Chebyshev-tau technique.

### (a) Algorithm

Solution of the generalized eigenvalue problem (3.2) yields a discrete, finite spectrum of eigenvalues {*c*_{0}, *c*_{1}, *c*_{2}, …} and eigenfunctions {{*u*_{0}(*x*), *v*_{0}(*x*)}, {*u*_{1}(*x*), *v*_{1}(*x*)}, …} where *c*_{j} corresponds to the pair {*u*_{j}(*x*), *v*_{j}(*x*)} and we use the ordering . Referring to the role of *c* in the form (2.35), we find that instability occurs when , and so the critical disturbance to track is modelled by {*u*_{0}(*x*), *v*_{0}(*x*)}.

It is clear from equations (2.30) and (2.31) that the entries of matrices *A* and *B* depend on our wavenumber *a*, Reynolds number *Re*, our position in the parameter space *Γ*={*δ*, *W*, *f*, *λ*_{1}, *λ*_{δ}} and our choice of *n* (see (2.35) for the role of *n*). With our parameters chosen, we iterate over *a* and *Re* tracking the ‘neutral curves’ along which , which separate regions of linear stability from those of instability. We do this first with *n*=0, and then increment *n* until we are at the computational limits of our method. As well as finding the neutral curves, we are particularly interested in the value , i.e. the minimum Reynolds number guaranteeing the onset of instability.

There is the need to move about *Γ* in a methodical way. The matter is simplified somewhat by choosing only to consider the case *λ*_{1}=*λ*_{δ}=*λ*, but we further simplify the problem by investigating the effects of each of *W*, *f*, *λ* in isolation from each other, in each case for at least the two values *δ*=0.4 (i.e. a thin thread) and *δ*=0.7 (a large diameter thread). The reasons for these specific values of *δ* will be made clear in our results.

## 4. Results

Figure 1 shows our results for how *Re*_{crit} varies with *δ*, for *n*=0, 1, …, 5, when *W*=*f*=*λ*=0. The *n*=0 curve agrees with the results of Walton (2004), whose non-dimensionalization we share, and we remind the reader that this is not the non-dimensionalization employed by Mott & Joseph (1968), who found a monotonic relation between *δ* and the critical Reynolds number based on their definition. However, we also include curves corresponding to the asymmetric cases of *n*>0. The behaviour of *Re*_{crit} at *δ*=0 and 1 is not clear, owing to the method of becoming numerically unstable at high Reynolds numbers. Interestingly, we find that axisymmetric perturbations appear to dominate only as *δ* approaches 1, i.e. as we tend towards plane Poiseuille flow, while the *n*=4 and 5 perturbations never dominate (table 1).

For all aspect ratios in the range 0.1<*δ*<0.8, we found the critical Reynolds number to be the minimum point of a unique neutral curve, qualitatively similar to that for plane Poiseuille flow, see the illustrations in Drazin & Reid. As *a* is increased, the wavenumber *a* at which the critical Reynolds number occurs also increases. No neutral curves could be found for *n*≥6.

### (a) Stabilizing effect of slip boundary conditions

Figure 1 shows the region 0.4<*δ*<0.7 to be the most interesting interval in terms of the behaviour of *Re*_{crit}. Not only does the periodicity of the least stable mode increase from *n*=1 to 3, but the ordering of the second and third least stable modes by periodicity is also seen to change. Therefore, in investigating the effect of our parameter space *Γ*, it is reasonable to restrict our computations to the cases *δ*=0.4 and 0.7, and values in between.

The linear analyses of Spille *et al*. (2000), Lauga & Cossu (2005), and the linear and nonlinear analyses of Webber & Straughan (2006) have shown boundary slip to have a stabilizing effect on Poiseuille flow in the plane. That is, *Re*_{crit} becomes higher with increasing slip length *λ*. This was in contrast to convective instabilities where increasing *λ* was found to have a destabilizing effect (see Webber 2006). Our results for the thread-annular case support the assertion that boundary slip is stabilizing in the case of parallel flows.

Having set *W*=*f*=0, figure 2 shows the critical Reynolds number to be a monotone increasing function of *λ* for each value of *n*. *λ* is also seen to influence the ordering of the second and third least stable perturbations: in the *δ*=0.4 case, the second least stable disturbance is seen to switch from being the *n*=2 mode to the axisymmetric mode at *λ*∼0.35. A similar reordering is seen when *δ*=0.6 around *λ*∼0.38.

Once more, we find that the neutral curve is unique and resembles the classic neutral curve associated with Poiseuille flow in the plane, see for example figure 3. As we increase *λ*, the critical wavenumber increases, and the wavenumber at which *Re*_{crit} occurs decreases. However, in figure 3 we find that the neutral curve begins to close by *λ*=0.011. Though we did not observe this for other values of *n* and *δ*, it implies that at finite *λ* the system becomes linearly stable to at least some types of perturbation.

### (b) Stabilizing effect of thread velocity *W*

We investigate the effect of *W*, setting *f*=*λ*=0. The most important work on this problem is that of Walton (2005), who performed calculations for the two cases *δ*=0.4 and 0.55, for disturbances with periodicities *n*=0, 1, 2. In obtaining our results on the effect of varying *δ* and *λ* (see above), we can corroborate Walton's assertion that in the case *W*=0 the neutral curve is always unique.

For the axisymmetric case (*n*=0) and *δ*=0.55, Walton found that for *W*>0 a stable intrusion appeared in the neutral curve, which grew with increasing *W* until the neutral curve was split into two separate curves. The ‘upper’ curve (in the domain of higher wavenumbers) closed up at finite *W*, while the lower curve retreated to infinity, such that ultimately this specific system became linearly stable to axisymmetric perturbations. He did not find this behaviour when setting *δ*=0.4. Instead, he reports of a unique neutral curve that closes with increasing *W*.

However, in our computations for *δ*=0.4 and *n*=0, we find that a stable intrusion does indeed occur for *W*∼*O*(10^{−3}), as seen in figure 4*a*–*d*. The smaller, closed curve that results from this stable intrusion vanishes for *W*>0.00799, and as we increase *W* further the remaining curve, which was identified by Walton (2005), closes (see figure 4*e*,*f*).

We observed a similar behaviour for *δ*=0.7—the neutral curve is unique when *W*=0; when *W*>0, a stable intrusion splits the neutral curve from right to left (in the orientation of figure 4); for large enough *W*, the neutral curve is completely split into two separate neutral curves. However, when *δ*=0.7 the smaller curve, observed to close up at finite *W*, was ‘above’, i.e. it was to be found at higher wavenumbers, while the larger ‘lower’ curve persisted for high values of *W*.

Walton (2005) also performed computations for the cases *δ*=0.4, *n*=1, 2. He reported that for both types of asymmetric perturbation, the neutral curve was observed to be unique. However, the results of our computations, see figure 5, again show the formation of a stable intrusion that splits the neutral curve. The smaller of the resulting curves is eventually seen to close up, but that larger curve tracked by Walton persists for large *W*.

The phenomenon of stable intrusion followed by curve splitting in the case of asymmetric perturbations was also observed for *δ*=0.7, with *n*=1. However, we did not see this when *δ*=0.4 and *n*=2, or when *δ*=0.7 and *n*≥2. Reasons for this include the neutral curve being very thin, making searching for neutral stability by iterating in the *a* and *Re* directions problematic and the numerical instability found at higher Reynolds numbers and for high values of *n*.

When curve splitting did occur, the critical Reynolds number was always found to lie on the larger of the two curves. Our numerical results would therefore indicate that when a smaller neutral curve exists, it is at higher Reynolds numbers (in the computationally stiff region), and so they do not affect the systems' transition to instability. Therefore, we present the results from tracking what we believe to be the critical Reynolds number in figure 6.

The results show that for *W* close to zero, increasing the thread velocity has a stabilizing effect. As *W* continues to increase, *Re*_{crit} begins to fall slightly, before increasingly rapidly, such that the overall influence of *W* is to stabilize the system. Interestingly, this behaviour is seen in the curves in both axisymmetric and asymmetric modes. For *δ*=0.4 the critical perturbation is consistently the *n*=1 mode, whereas when *δ*=0.7 the critical mode changes with *W*, see table 2.

### (c) Destabilizing effect of rotation

Finally, we investigate the effect of allowing the thread to rotate. We remind the reader of the role of *f* in the boundary conditions in (2.8). *f* is the angular frequency with which our thread of radius *δ* spins around the *z*-axis, such that points on the threads surface have velocity (when *W*=0), and we showed that *f*>0 introduced a new component into the steady state.

For both aspect ratios, our results for axisymmetric perturbations showed the neutral curve to be unique. Increasing *f* had the effect of slowly decreasing the value of *Re*_{crit}, so that rotation was seen to be destabilizing, see figure 7. However, this was not true in the case of asymmetric perturbations.

For an aspect ratio of *δ*=0.7, and asymmetric perturbations *n*=1, 2, 3, increasing *f* caused a stable intrusion to appear, splitting the neutral curve into two separate curves. The minimum Reynolds number of these two was tracked and shown to increase with *f*, stabilizing the flow. However, when *f*>0 a new, separate neutral curve appears in the region of very small wavenumbers. As *f* increases this neutral curve creeps from right to left, until eventually it contains the critical Reynolds number, see figure 8. From this point on, *Re*_{crit} is a monotone decreasing function of *f*. This behaviour was also seen when *δ*=0.4 for the *n*=1 perturbations.

For other combinations of aspect ratio and *n*>0, we observed this new, destabilizing neutral curve appear, but we did not find curve splitting by stable intrusion. This is because of the following: often the neutral curve becomes very thin, and so curve splitting can easily be missed by our iteration over *a* and *Re*; the method becomes unstable at high Reynolds numbers, especially as *f* is increased; the new neutral curve expands and obscures the other neutral curve before curve splitting takes place (figure 9).

The eventual dominance of this new neutral curve is interesting, in that it implies a sudden change in the system's instability. At a critical value of *f*, long-wave (i.e. small wavenumber) perturbations suddenly become the most unstable form of disturbance.

## 5. Conclusions and further work

Our results on the effects of *λ* and *f* represent entirely new and original work. In this problem we have found yet another parallel flow problem for which boundary slip is stabilizing, and the destabilizing influence of thread rotation may be of direct importance in the use of thread injection.

Although treatment of a constant thread speed *W* is already present in the literature, we have been thorough in our computations and spotted behaviour missed by previous authors. For the first time, the critical Reynolds number is shown as a function of *W*, for both asymmetric and axisymmetric perturbations.

One of the main limitations of our results is the computer RAM available to us, which directly impacts upon the number of polynomials we can use in our numerical method. Equation (3.2) was found to be very stiff, and even when using the maximum number of Chebyshev polynomials possible, the method's numerical instability at high Reynolds numbers, or for high values of *n*, or large *f*, was problematic.

There is much scope for future work on this subject. For instance, with regards to the model: allowing the thread to flex or not be concentric with the outer cylinder (as suggested by Walton 2005); allowing the thread itself to be a porous medium; consideration of the thread's finite length in the interior of the hypodermic needle.

## Acknowledgments

This work was supported by a research studentship of the Engineering and Physical Sciences Research Council. The author would like to thank Prof. B. Straughan for his helpful advice.

## Footnotes

- Received April 10, 2007.
- Accepted November 19, 2007.

- © 2008 The Royal Society