## Abstract

The surgical technique of thread injection of medical implants is modelled by the axial pressure-gradient-driven flow between concentric cylinders with a moving core. The linear stability of the flow to both axisymmetric and asymmetric perturbations is analysed asymptotically at large Reynolds number, and computationally at finite Reynolds number. The existence of multiple regions of instability is predicted and their dependence upon radius ratio and thread velocity is determined. A discrepancy in critical Reynolds numbers and cut-off velocity is found to exist between experimental results and the predictions of the linear theory. In order to account for this discrepancy, the high Reynolds number, nonlinear stability properties of the flow are analysed and a nonlinear, equilibrium critical layer structure is found, which leads to an enhanced correction to the basic flow. The predictions of the nonlinear theory are found to be in good agreement with the experimental data.

## 1. Introduction and formulation of the governing equations

### (a) Thread injection

Thread injection is a recently devised technique enabling porous medical implants to be placed inside the body in a minimally invasive manner, thus reducing surgical trauma. The thread is stored on a spool and injected within a fluid by applying an axial pressure gradient, with the thread velocity controlled by a motor (figure 1). In a recent experimental paper, Frei *et al*. (2000) modelled the thread injection process by using a cylindrical rubber filament to represent the thread and a steel outer pipe filled with water. On comparing their results with the exact Navier–Stokes solution for sliding flow through an annulus (equation (1.1) below), they discovered that the force on the thread measured in the experiments was always significantly less than theory predicts. This observation provides the impetus for the current study in which we investigate whether this discrepancy could be caused by a linear or nonlinear instability of the basic flow.

The linear stability of this flow has been studied computationally by Sadeghi & Higgins (1991) and Gittler (1993). There is some disagreement between their results, particularly concerning the number of neutral curves that coexist at given values of thread radius and velocity. In order to resolve these differences, we begin in §2*a* by analysing the linear stability equations asymptotically at large Reynolds number *R*. We then present some finite *R* computations for both axisymmetric and non-symmetric disturbances. In §3, we discuss the nonlinear stability of the flow at large *R* and give predictions for the force on the thread based on this theory. Possible future work is discussed in §4.

### (b) The governing equations, non-dimensionalization and basic flow

The governing equations are the non-dimensional, incompressible Navier–Stokes equations, with Reynolds number *R* based on the constant axial pressure gradient applied to the flow and the fixed radius of the outer cylinder. Our basic state is the axial flow between concentric cylinders with *r*=1 representing the outer surface and *r*=*δ* the inner surface. The inner cylinder is moving in the axial direction with non-dimensional velocity *V*. We seek a steady, unidirectional solution * u*=(

*U*

_{0}(

*r*),0,0) to the Navier–Stokes equations in cylindrical polar coordinates (

*x*,

*r*,

*θ*) and find(1.1)

In this paper, we study the linear and nonlinear stability of this basic solution for various values of *V* and *δ*. We will concentrate on *V*>0, as this is the regime of interest for medical applications. The case *V*<0 is dealt with in Walton (2004).

## 2. Linear stability

A small travelling wave disturbance of non-dimensional amplitude *Δ* is superimposed upon the basic flow and is expressed as the real part of(2.1)

A temporal approach to the stability problem is taken, allowing the wavespeed *c* to be complex, while the axial wavenumber *α* is required to be real. The governing equations for the disturbance are obtained by adding the perturbation to the basic flow, substituting into the Navier–Stokes equations and equating coefficients of *Δ*, leading to a set of linear equations. It is easy to show, using a form of Rayleigh's inflection point theorem, that in the inviscid limit *R*=∞, the flow is linearly stable. Thus, any neutral stability curves would be expected to close up, either at finite Reynolds number or as *R*→∞.

In the situation where the perturbations are axisymmetric in nature (*N*=0), the disturbance equations can be manipulated into a fourth-order Orr–Sommerfeld type problem:(2.2)where , with the boundary conditions *ϕ*=*ϕ*′=0 on *r*=*δ*, *r*=1. The solution of (2.2) for given values of *V* and *δ* is a numerical one in general, but some analytical progress can be made by assuming that the Reynolds number is large. This asymptotic analysis will prove useful when we seek to explain the behaviour of the solutions at finite values of *R* (§2*b*).

### (a) Asymptotic solutions for axisymmetric disturbances

When the thread velocity *V*=0 there are two distinguished scalings, each yielding one solution, that describe the behaviour on the lower and upper branches of the linear neutral stability curve at sufficiently high Reynolds number. These scalings remain intact for a range of non-zero *V*, but extra solutions are created, suggesting that for certain flow configurations, multiple neutral curves exist. This is confirmed by the finite *R* numerical calculations we present in §2*b*. First, we discuss the asymptotic solutions and their behaviour as *V* is increased. The analysis presented below is similar to Cowley & Smith's (1985) study of plane Poiseuille–Couette flow, and can be found in more detail in Walton (2004).

#### (i) Lower branch modes

The lower branch structure, consisting of two viscous wall layers and an inviscid core, remains essentially unchanged up to and including values of *V* of O(*R*^{−2/7}). At this stage the appropriate scalings are(2.3)

Since, in practice, we are interested in the stability of this flow at O(1) values of *V*, the limit *V*_{0}→∞ is relevant here. In this limit, it can be shown that for *δ* in the range 0.92<*δ*<1, there are four solutions of the lower branch eigenrelation for which *c*_{0} remains O(1). The number of solutions reduces to two for radius ratios 0.52<*δ*<0.92 and there are no solutions of this form for *δ*<0.52. The wavelength of each of these modes increases with increasing *V*_{0}, with as *V*_{0}→∞ (from asymptotic analysis of the relevant eigenrelation). When *V* is increased to O(1), the wavelength of these modes has therefore increased to O(*R*) (in view of the scaling (2.3)), and their further development must be calculated numerically (see §2*a*(iv) below). In addition to these modes, for all values of *δ*, there is a mode whose wavelength shortens, with as *V*_{0}→∞, while the corresponding wavespeed approaches the thread velocity from below. For this mode, the critical layer therefore moves away from the upper wall as *V* increases, while the lower critical layer remains embedded within its viscous layer. Consequently, the mode becomes a mixture of an upper branch and lower branch mode: we refer to it as a ‘hybrid’ mode and discuss it further in §2*a*(iii) below.

#### (ii) Upper branch modes

The upper branch structure, in which the critical and viscous layers are separated, remains valid up to and including thread velocities of O(*R*^{−2/11}). At that stage the appropriate scalings are

It can be shown that a unique solution exists for a given flow configuration. Typically, *α*_{0} increases with increasing *δ*, and as *V*_{0} is increased from zero, the value of *α*_{0} rises slightly, falls and then increases monotonically with , *c*_{0}→*V*_{0} as *V*_{0}→∞. Once again, we see that the wavespeed is approaching the thread velocity, indicating that from the upper branch perspective, the lower critical layer is moving towards the wall. Again, we expect that with further increase in *V* this mode will evolve into the hybrid type mentioned above.

#### (iii) Hybrid modes

By examining the breakdown of the lower and upper branch solutions as *V* increases, it is possible to derive the appropriate scalings for the hybrid mode. These are, to leading order:

It is straightforward to demonstrate that in the limit , the solutions of the hybrid eigenrelation match back to the lower and upper branch modes discussed earlier. Up to a critical value of , there are typically two solutions to this eigenrelation, with no solutions beyond this value. The neutral curve corresponding to the hybrid mode therefore experiences ‘cut-off’ at a value of *V* of O(*R*^{−2/13}) when *R* is large. Although this suggests that such modes will not be present at O(1) values of *V*, our numerical calculations (§2*b*) show that at finite *R* this mode persists to larger values of *V* than might have been expected from the asymptotic analysis.

#### (iv) Long wave modes

We have seen that as *V* is increased, the upper branch solution and one of the lower branch modes evolve onto the hybrid scaling and disappear completely on that scaling as *V* is further increased. The ultimate fate of the other four solutions generated on the lower branch scaling remains to be addressed. As remarked in §2*a*(i) above, when *V* becomes O(1), these modes acquire long waves properties with *α*∼O(*R*^{−1}). Their subsequent evolution as *V* is increased further is determined by taking the limit *R*→∞, but allowing *αR* to remain O(1). Thus, the modes are governed by (2.2) with terms proportional to *α*^{2} and *α*^{4} set to zero; we will refer to this as the ‘reduced Orr–Sommerfeld equation’. Numerical solutions of this equation show that one pair of modes disappears at very small *V* on this scale with the corresponding value of *αR* being large, suggesting that these modes only exist at very large *R* and have no physical relevance. The second set of modes experiences cut-off at much larger values of *V* (and smaller *αR*). For each value of *δ*, there is a corresponding *V*_{c}, but below *δ*=0.428, there are no modes at all on this scaling. However, this does not mean that the flow is completely linearly stable below this radius ratio as we shall see below when we discuss the finite Reynolds number calculations.

### (b) Solutions at finite Reynolds number

We now discuss the results of a numerical solution of the linear disturbance equations. A Chebyshev collocation approach of the form described in Schmid & Henningson (2001) for *N*=0 and Khorrami *et al*. (1989) for *N*≠0 was used to approximate the continuous system and an arclength continuation approach suggested by Dr G. Moore (2003, personal communication) was used to follow the trajectory of each neutral curve. First we discuss the axisymmetric results.

#### (i) Axisymmetric calculations

Here, we concentrate on two representative radius ratios: *δ*=0.55 and *δ*=0.4. For the larger of these values of *δ*, the behaviour of the neutral stability regions as *V* is increased is as follows (see figure 2). First, for *V*=0, there is a unique neutral curve (with lower and upper branches described by the asymptotic analysis of §2*a*(*i*),(*ii*)). At small non-zero *V*, a stable intrusion forms and spreads from right to left, eventually penetrating almost to the nose of the curve and slicing it into two, approximately equal sized, smaller curves. The top curve is described by the hybrid scaling discussed in §2*a*(*iii*), while the lower curve is governed by the reduced Orr–Sommerfeld equation discussed in §2*a*(*iv*). As *V* is increased further, the hybrid curve thins, closes up at finite *R*, and eventually shrinks to a point. By contrast, the lower curve persists over a larger range of *V*, eventually retreating to *R*=∞ at a finite value of *V*, which can be predicted by the analysis of §2*a*(*iv*). Beyond this value of *V*, the flow is completely linearly stable to axisymmetric disturbances. Neutral curves for the lower radius ratio *δ*=0.4 are shown in figure 3. In this case, no stable intrusion forms and the unique curve for *V*=0 simply evolves onto the hybrid scaling, eventually closing up and disappearing in a similar manner to that described above.

#### (ii) Calculations for non-symmetric modes

In this case, no asymptotic analysis for *N*∼O(1) seems to be possible, and the numerical results suggest that there is always a unique neutral curve for a particular flow configuration. Neutral curves are shown in figure 4 for the case *δ*=0.4. In figure 4*a*, we examine the stability of the *N*=1 mode, where we see that the wavelength of the instability increases rapidly as *V* is increased. The *V*=0 solution is most unstable with a critical Reynolds number slightly smaller and a cut-off velocity significantly larger than the axisymmetric mode. Figure 4*b* shows the stability of the *N*=2 mode as *V* is increased. The critical Reynolds number *R*_{crit} is slightly larger than for *N*=1 and the cut-off velocity is much smaller. Additional numerical studies indicate that the flow is least stable when *δ*≃0.2, and for a given *δ*, increasing *V* always leads to an increase in *R*_{crit}. The cut-off velocities observed here can be calculated accurately (as in the *N*=0 case) by considering the limit *R*→∞, *αR*→O(1) of the disturbance equations and below a thread radius of about 0.5, it appears that this limit describes the instability even when *V*=0. For non-symmetric modes, there appears to be no analogue of the hybrid mode and all of the neutral curves retreat to infinity as *V* is increased, rather than closing up at finite Reynolds number.

### (c) Summary of the linear results

Using a combination of asymptotic and numerical analyses, we have investigated the linear stability of thread-annular flow to axisymmetric and asymmetric disturbances. For *N*=0, multiple regions of instability coexist and are governed by different scalings. For non-symmetric disturbances, there is a unique neutral curve for each *N*, with the *N*=1 mode being the most unstable. Every mode has a cut-off velocity *V*_{c}, above which the flow is linearly stable at all *R*, and the *V*=0 configuration is always the most unstable for a given thread radius *δ*. For the *N*=1 mode, the flow achieves its lowest critical Reynolds number (approximately 4×10^{5}) when *δ*≃0.2. If we compare the figures for *R*_{crit} with those obtained in the experiments of Frei *et al*. (2000), we see that in terms of our non-dimensionalization, they reported instability at *R*≃3×10^{4}. Therefore, there is an order of magnitude disagreement between linear theory and experiment, suggesting that nonlinear effects play a key role in the transition process, presumably in a similar way to that for the flow through a single pipe. We investigate this further in §3.

## 3. Nonlinear stability at large Reynolds number

We investigate the nonlinear stability of thread-annular flow at large Reynolds number by seeking a nonlinear travelling-wave type solution of the form found by Smith & Bodonyi (1982*a*) for the flow through a single pipe. A more detailed version of the following analysis can be found in Walton (2003).

### (a) Determination of the disturbance amplitude

As in the linear analysis, we superimpose a disturbance of the form (2.1), but now we concentrate on a particular disturbance size, *Δ*=*R*^{−1/3}*A*, with the O(1) value of *A* to be determined as a function of *α*, *V* and *δ*. The scaling for *Δ* proposed here is implied by previous studies of nonlinear equilibrium critical layers (e.g. Smith & Bodonyi 1982*a*,*b*), and ensures that the small phase shifts induced across the critical and wall layers are of the same order of magnitude. This particular choice for *Δ* also leads to a relatively large mean-flow correction to the basic streamwise flow of the form *R*^{−1/6}*u*_{M}(*r*), with a similar correction to the azimuthal velocity. The flow in the majority of the pipe is inviscid, with the pressure perturbation *P* satisfying Rayleigh's equation subject to Neumann conditions at the walls, and the condition of an asymptotically small phase shift at the critical level *r*_{c} (where *U*_{0}(*r*_{c})=*c*), so that terms of the form (ln(*r*_{c}−*r*)exp(i*ξ*)) for *r*<*r*_{c} are replaced by ((ln(*r*−*r*_{c})+i*R*^{−1/2}*Φ*_{CL})exp(i*ξ*)) for *r*>*r*_{c}, where the O(1) quantity *Φ*_{CL} is to be found in terms of *A* from the critical layer analysis. The Rayleigh problem can be solved by a Runge–Kutta method coupled with an iterative approach in which for given *V*, *δ* the value of *r*_{c} is at first guessed and then updated using Newton's method. Once the solution for *r*_{c} is obtained, the values of *P* at *r*=*δ*, 1 can be found, as can the leading order wavespeed *c*. Turning now to the determination of the phase shift, after an extensive analysis of the critical layer, we find a phase shift-amplitude relation of the form *Φ*_{CL}=*A*^{−3/2}*f*(*α*, *N*, *δ*, *V*). This phase shift must be in tune with that induced across the viscous wall layer. Setting the phase shifts equal, we obtain the neutral wave amplitude as follows:(3.1)Here, *C*^{(1)} and *C*^{(2)} are numerical constants arising from the finite parts of certain integrals with *C*^{(1)}≃−5.516, *C*^{(2)}≃−0.1564, while *P*(1), *P*(*δ*) are known from the solution of the Rayleigh problem in the core. In figure 5, we show the amplitude of the *N*=1 mode versus *δ* for various values of injection velocity *V* for the particular case *α*=1. We also show the corresponding dependence of *r*_{c} upon *δ*. Similar plots are obtained for other values of *α*. We can see from figure 5*a* that the nonlinear cut-off velocity for *δ*=0.4 is *V*_{c}≃0.3, considerably in excess of the linear cut-off found earlier. For larger values of *N*, the behaviour is qualitatively similar but the amplitudes involved are much smaller, suggesting that the *N*=1 mode is the most dangerous in practice, as was found to be the case for linear disturbances.

### (b) The mean-flow distortion and the force on the thread

In order to calculate the force on the thread, we need to find the correction to the mean flow induced by this nonlinear disturbance. From substitution into the Navier–Stokes equations, we find that the distortion is governed by . Although there is no wave forcing on the right-hand side of this equation, the disturbance makes its presence felt via jump conditions on *u*_{M}, across the critical layer. The critical layer analysis (Walton 2003) indicates that the appropriate solution for *u*_{M} is(3.2)where the quantities *M*_{1} and *M*_{2} can be written down explicitly in terms of the parameters *A*, *δ*, *V*, *α*, *N* and *r*_{c}. The force on a thread of length *L* owing to viscous stresses and the axial pressure gradient is given in non-dimensional terms by *F*/4*RL*=*πδu*_{r}(*δ*)/2+*πδ*^{2}. Taking *u*=*U*_{0}+*R*^{−1/6}*u*_{M} and using (3.2) and (3.1), we obtain(3.3)

In figure 6, we plot the resistive force calculated from equation (3.3) (with *N*=1) as a function of *δ* for two different injection velocities *V*=0.08 and *V*=0.2. These correspond to values used in the experiments of Frei *et al*. (2000). The experiments were performed at a flux based Reynolds number Re=100, where Re is related to *R* via Re=2*RQ*/*π*(1+*δ*), with *Q*, the dimensionless flux of fluid through the annulus. It can be seen that our results are relatively insensitive to changes in the axial wavenumber *α* and that, in all cases, the thread force is reduced by the presence of the instability. The agreement with the experimental data is encouraging, indicating that a nonlinear instability of this form could indeed be responsible for the discrepancy between theory and experiment.

## 4. Discussion and future work

We have investigated the linear stability of thread-annular flow to axi- and non-symmetric disturbances both asymptotically and at finite Reynolds numbers. The results show that non-symmetric disturbances are more dangerous; the most unstable configuration is when the ratio of inner to outer cylinder is approximately 0.2, and that increasing the thread velocity raises the critical Reynolds number, rendering the flow more stable. For all configurations, there is a cut-off velocity beyond which the flow is stable to all infinitesimal disturbances. These cut-off velocities are generally less than those used in the experiments of Frei *et al*. (2000) where instability is reported. The critical value of Reynolds number in the experiments is also often found to be an order of magnitude lower than that predicted by linear theory. Thus, it appears that the discrepancy between theory and experiment cannot be the result of a linear instability. We then considered the nonlinear instability using a high Reynolds number approach in which the perturbation takes the form of a single neutral mode across the majority of the annulus. The instability is predominantly governed by inviscid dynamics, but there exists a strongly nonlinear critical layer where viscous effects determine the phase shift and hence the disturbance amplitude. This structure leads to an enhanced mean-flow correction and therefore makes a significant contribution to the force on the thread. The results of this theory are compared with the experiments and encouraging agreement is found.

It should be noted that the simple model of axial flow between concentric cylinders is a great simplification of the actual thread injection problem. In practice, there are many other effects that need to be taken into account and thus our results should be regarded as merely a first step in describing the experimentally observed phenomena. In reality, the thread will possess some flexibility and its position within the outer cylinder may not always be concentric. Both of these factors could affect the instability process considerably. There is also the issue of whether the thread is sufficiently long for the assumption of a spatially fully developed basic flow to be valid, and whether the injection process occurs over a long enough period of time for temporal variations in *U*_{0} to be neglected. It is hoped that further studies in this area can incorporate these extra effects.

## Footnotes

One contribution of 19 to a Theme ‘New developments and applications in rapid fluid flows’.

- © 2005 The Royal Society