# [

## Abstract

Motivated by complex multi-fluid geometries currently being explored in fibre-device manufacturing, we study capillary instabilities in concentric cylindrical flows of fluids with arbitrary viscosities, thicknesses, densities, and surface tensions in both the Stokes regime and for the full Navier–Stokes problem. Generalizing previous work by Tomotika (), Stone & Brenner (, equal viscosities) and others, we present a full linear stability analysis of the growth modes and rates, reducing the system to a linear generalized eigenproblem in the Stokes case. Furthermore, we demonstrate by Plateau-style geometrical arguments that only axisymmetric instabilities need be considered. We show that the case is already sufficient to obtain several interesting phenomena: limiting cases of thin shells or low shell viscosity that reduce to problems, and a system with competing breakup processes at very different length scales. The latter is demonstrated with full 3-dimensional Stokes-flow simulations. Many cases remain to be explored, and as a first step we discuss two illustrative cases, an alternating-layer structure and a geometry with a continuously varying viscosity.

Stability analysis of concentric shells]Linear stability analysis of capillary instabilities for concentric cylindrical shells

X. Liang, D. S. Deng, J.-C. Nave, and S. G. Johnson]X.\nsL\lsI\lsA\lsN\lsG,\nsD.\nsS.\nsD\lsE\lsN\lsG,J.\ls-\lsC.\nsN\lsA\lsV\lsE\ns and S\lsT\lsE\lsV\lsE\lsN\nsG.\nsJ\lsO\lsH\lsN\lsS\lsO\lsN

apillary instability, cylindrical shells

## 1 Introduction

In this paper, we generalize previous linear stability analyses Lord Rayleigh (1879, 1892); Tomotika (1935); Stone & Brenner (1996); Chauhan et al. (2000) of Plateau–Rayleigh (capillary) instabilities in fluid cylinders to handle any number () of concentric cylindrical fluid shells with arbitrary thicknesses, viscosities, densities, and surface tensions. This analysis is motivated by the fact that experimental work is currently studying increasingly complicated fluid systems for device-fabrication applications, such as drawing of microstructured optical fibres with concentric shells of different glasses/polymers Hart et al. (2002); Kuriki et al. (2004); Pone et al. (2006); Abouraddy et al. (2007); Sorin et al. (2007); Deng et al. (2008) or generating double emulsions Utada et al. (2005); Shah et al. (2008). Although real experimental geometries may not be exactly concentric, we show that surface tension alone, in the absence of other forces, will tend to eliminate small deviations from concentricity. We show that our solution reduces to known results in several limiting cases, and we also validate it with full 3-dimensional Stokes-flow simulations. In addition, we show results for a number of situations that have not been previously studied. For the limiting case of a thin shell, we show a connection to the classic single-cylinder and flat-plane results, consistent with a similar result for air-clad two-fluid jets Chauhan et al. (2000). In a three-fluid system, we exhibit an interesting case in which two growth modes at different wavelengths have the same effective growth rate, leading to competing breakup processes that we demonstrate with full 3-dimensional Stokes-flow simulations. We also consider some many-layer cases, including a limiting situation of a continuously varying viscosity. Using a simple geometrical argument, we generalize previous results Plateau (1873); Lord Rayleigh (1879); Chandrasekhar (1961) to show that only axial (not azimuthal) instabilities need be considered for cylindrical shells. Numerically, we show that the stability analysis in the Stokes regime can be reduced to a generalized eigenproblem whose solutions are the growth modes, which is easily tractable even for large numbers of layers. Like several previous authors Tomotika (1935); Stone & Brenner (1996); Gunawan et al. (2002, 2004), we begin by considering the Stokes (low-Reynolds) regime, which is consistent with the high viscosities found in drawn-fibre devices Abouraddy et al. (2007); Deng et al. (2008). In Appendix C, we generalize the analysis to the full incompressible Navier-Stokes problem, which turns out to be a relatively minor modification once the Stokes problem is understood, although it has the complication of yielding an unavoidably nonlinear eigenproblem for the growth modes. Semi-analytical methods are a crucial complement to large-scale Navier-Stokes simulations (or experiments) in studying capillary instabilities, since the former allow rapid exploration of wide parameter regimes (e.g. for materials design) as well as rigorous asymptotic results, while the latter can capture the culmination of the breakup process as it grows beyond the linear regime.

Capillary instability of liquid jets has been widely studied Lin (2003); Eggers & Villermaux (2008) since Plateau (1873) first showed that whenever a cylindrical jet’s length exceeds its circumference, it is always unstable due to capillary forces (surface tension). Plateau used simple geometrical arguments based on comparing surface energies before and after small perturbations. Lord Rayleigh introduced the powerful tool of linear stability analysis and reconsidered inviscid water jets Lord Rayleigh (1879) and viscous liquid jets Lord Rayleigh (1892). In linear stability analysis, one expands the radius as a function of axial coordinate in the form , where , is a wavelength of the instability, and is an exponential growth rate. Given a geometry, one solves for the dispersion relation(s) and considers the most unstable growth mode with the growth rate to determine the time scale of the breakup process. The wavelength corresponding to has been experimentally verified to match the disintegration size of liquid jets Eggers & Villermaux (2008). By considering the effect of the surrounding fluid, Tomotika (1935) generalized this analysis to a cylindrical viscous liquid surrounded by another viscous fluid, obtaining a determinant equation for the dispersion relation. (Rayleigh’s results are obtained in the limit of vanishing outer viscosity.) A few more limiting solutions have been obtained Meister & Scheele (1967); Kinoshita et al. (1994) by generalizing Tomotika’s approach to non-Stokes regimes. Beyond the regime of a single cylindrical jet, Stone & Brenner (1996) analysed the three-fluid () Stokes cylinder problem, but only for equal viscosities. Chauhan et al. (2000) analysed the case where the inner two fluids have arbitrary viscosities and the outermost fluid is inviscid gas, taking into account the full Navier–Stokes equations. Gunawan et al. (2002, 2004) considered an array of identical viscous cylinders (in a single row or in a triangular configuration).

Much more complicated, multi-fluid geometries are now being considered in experimental fabrication of various devices by a fibre-drawing process Hart et al. (2002); Kuriki et al. (2004); Pone et al. (2006); Abouraddy et al. (2007); Sorin et al. (2007); Deng et al. (2008). In fibre drawing, a scale model (preform) of the desired device is heated to a viscous state and then pulled (drawn) to yield a long fibre with (ordinarily) identical cross-section but much smaller diameter. For example, concentric layers of different polymers and glasses can be drawn into a long fibre with submicron-scale layers that act as optical devices for wavelengths on the same scale as the layer thicknesses Abouraddy et al. (2007). Other devices, such as photodetectors Sorin et al. (2007), semiconductor filaments Deng et al. (2008, 2010), and piezoelectric pressure sensors Egusa et al. (2010) have similarly been incorporated into microstructured fibre devices. That work motivates greater theoretical investigation of multi-fluid geometries, and in particular the stability (or instability time scale) of different geometries is critical in order to predict whether they can be fabricated successfully. For example, an interesting azimuthal breakup process was observed experimentally Deng et al. (2008) and has yet to be explained ?; however, we show in this paper that azimuthal instability does not arise in purely cylindrical geometries and must stem from the rapid taper of the fibre radius from centimetres to millimetres (the drawn-down “neck”), or some other physical influence. [For example, there may be elastic effects (since fibres are drawn under tension), thermal gradients at longer length scales in the fibre (although breakup occurs at small scales where temperatures are nearly uniform), and long-range (e.g., van der Waals) interactions in submicron-scale films.] Aside from fibre drawing, recent authors have investigated multi-fluid microcapillary devices, in which the instabilities are exploited to generate double emulsions, i.e. droplets within droplets Utada et al. (2005). Because the available theory was limited to equal viscosities, the experimental researchers chose only fluids in that regime, whereas our paper opens the possibilities of predictions for unequal viscosity and more than three fluids.

We now formulate the mathematical problem that we solve, as depicted schematically in figure 1. The total number of viscous fluids is and the viscosity of the -th () fluid is . The surface-tension coefficient between the -th and fluid is denoted by . For the unperturbed steady state (figure b), we assume that the -th fluid is in a cylindrical shell geometry with outer radius and inner radius . The first () fluid is the innermost core and the -th fluid is the outermost one (extending to infinity), so we set and . To begin with, this system is analysed in the Stokes regime (low Reynolds number) and we also neglect gravity [in the large Froude number limit, valid for fibre-drawing Deng et al. (2011)], so the fluid densities are irrelevant. In Appendix C, we extend this analysis to the Navier–Stokes regime, including an inertia term for each layer (with density ). As noted above, linear stability analysis consists of perturbing each interface by a small sinusoidal amount , to lowest order in . Stokes’ equations are then solved in each layer in terms of Bessel functions, and matching boundary conditions yields a set of equations relating and . Although these equations can be cast in the form of a polynomial root-finding problem, similar to Tomotika (1935), such a formulation turns out to be ill-conditioned for large , and instead we formulate it in the Stokes regime as a generalized eigenproblem of the form , which is easily solved for the dispersion relations (with the corresponding eigenvectors yielding the relative amplitudes of each layer). In the Navier–Stokes regime, this becomes a nonlinear eigenproblem.

## 2 Azimuthal stability

For any coupled-fluids system of the type described in figure 1, a natural question to ask is whether that system is stable subject to a small perturbation. If an interface with area has surface energy , then a simple way to check stability is to compare surface energies (areas) for an initial configuration and a slightly deformed configuration with the same volume. In this way, it was shown that any azimuthal deformation is stable for a single cylindrical jet Lord Rayleigh (1879); Chandrasekhar (1961). Here, we employ a similar approach to demonstrate that the same property also holds for multiple concentric cylindrical shells. Note that this analysis only indicates whether a system is stable; in order to determine the time scale of an instability, we must use linear stability analysis as described in subsequent sections.

For the unperturbed system, we define the level-set function . defines the interface between the -th and -th fluids. Similarly, we define the level-set function for the perturbed interface (figure c) between the -th and -th fluids by

(1) |

Following the method of normal modes Drazin & Reid (2004), in the limit of small perturbations, a disturbed interface can be chosen in the form

(2) |

Assuming incompressible fluids in each layer so that volume is conserved (and assuming that the cylinder is much longer than its diameter so that any inflow/outflow at the end facets is negligible), we obtain a relation between and

(3) |

Let denote the surface area of in one wavelength . From equations (1) and (2), can be expressed in cylindrical coordinates as:

(4) |

Now we can compare the total interfacial energy between the unperturbed system and the perturbed system :

(5) |

From the surface-energy point of view, small perturbations will grow only if . Therefore, from (5), we can conclude that all the non-axisymmetric perturbations will be stable. There is one special case that needs additional consideration: if and in (5), the first term is zero, so one must consider the next-order term in order to show that this case is indeed stable (i.e. elliptical perturbations decay). Even more straightforwardly, however, corresponds to a two-dimensional problem, in which case it is well known that the minimal surface enclosing a given area is a circle.

## 3 Linear stability analysis

In the previous section, we showed that only axisymmetric perturbations can lead to instability of concentric cylinders. Now we will use linear stability analysis to find out how fast the axisymmetric perturbations grow and estimate the break up time scale for a coupled -layer system.

Here, we consider fluids in the low-Reynolds-number regime [valid for fibre-drawing Abouraddy et al. (2007); Deng et al. (2008)] and thus the governing equations of motion for each fluid are the Stokes equations Ockendon & Ockendon (1995). The full Navier–Stokes equations are considered in Appendix C. For the axisymmetric flow, the velocity profile of the -th fluid is , where is the radial component of the velocity and is the axial component of the velocity. The dynamic pressure in the -th fluid is denoted by . The Stokes equations Batchelor (1973) are

(6a) | |||

(6b) |

and the continuity equation (for incompressible fluids) is

(7) |

### 3.1 Steady state

Because of the no-slip boundary conditions of viscous fluids, without loss of generality, we can take the equilibrium state of the outermost fluid to be

(8) |

for , and of the -th fluid to be

(9) |

for . (In Appendix C, we generalize this to nonzero initial relative velocities for the case of inviscid fluids, where no-slip boundary conditions are not applied.)

### 3.2 Perturbed state

The perturbed interface, corresponding to the level set , with an axisymmetric perturbation, takes the form

(10) |

Similarly, the perturbed velocity and pressure are of the form

(11) |

Note that the Stokes equations (6) and continuity equation (7) imply that

(12) |

Substituting the third row of (11) into (12), we obtain an ordinary differential equation for

(13) |

Clearly, satisfies the modified Bessel equation of order zero in terms of . Therefore, we have

(14) |

where and are modified Bessel functions of the first and second kind ( is exponentially decreasing and singular at origin; is exponentially growing), and and are constants to be determined. Inserting into (6) and solving two inhomogeneous differential equations, we obtain the radial component of velocity

(15) |

and the axial component of velocity

(16) |

where and are constants to be determined. Imposing the conditions that the velocity and pressure must be finite at and , we immediately have

(17) |

### 3.3 Boundary conditions

In order to determine the unknown constants in each layer, we close the system by imposing boundary conditions at each interface. Let be the unit outward normal vector of interface and be the unit tangent vector. Formulae for and are given by equations (41) and (42) of Appendix A. First, the normal component of the velocity is continuous at the interface, since there is no mass transfer across the interface, and so

(18) |

For the at-rest steady state (8) and (9), this condition is equivalent (to first order) to the continuity of radial velocity:

(19) |

Second, the no-slip boundary condition implies that the tangential component of the velocity is continuous at the interface:

(20) |

(The generalization to inviscid fluids, where no-slip boundary conditions are not applied, is considered in Appendix C.) For the at-rest steady state (8) and (9), this is equivalent (to first order) to the continuity of axial velocity:

(21) |

Third, the tangential stress of the fluid is continuous at the interface. The stress tensor of the -th fluid in cylindrical coordinates for axisymmetric flow Kundu & Cohen (2007) can be expressed as

(22) |

The continuity of the tangential stress at the interface implies that

(23) |

The leading term of (23) leads to

(24) |

Fourth, the jump of the normal stress across the interface must be balanced by the surface-tension force per unit area. The equation for normal stress balance at the interface is

(25) |

where is the mean curvature of the interface. The curvature can be calculated directly from the unit outward normal vector of the interface by (see Appendix A). Substituting (22), (41), and (45) into (25), we have the following equation (accurate to first order in ) for the dynamic boundary condition:

(26) |

### 3.4 Dispersion relation

Substituting (11) and (14)–(16) into the boundary conditions [(19), (21), (24) and (26)] and keeping the leading terms, we obtain a linear system in terms of the unknown constants . After some algebraic manipulation, these equations can be put into a matrix form:

(27) |

where and are matrices given below.

(28) |

(29) |

Here, and denote the corresponding modified Bessel functions evaluated at . Combining the boundary conditions from all interfaces, we have the matrix equation

(30) |

for the undetermined constants , with

(31) |

where and are the second and fourth columns of and , and is the first and third columns of . To obtain a non-trivial solution of equation (30), the coefficient matrix must be singular, namely

(32) |

Since is zero except in its fourth row, only occurs in the th, th, , -th rows of . The Leibniz formula implies that equation (32) is a polynomial in with degree . Therefore, we could obtain the dispersion relation by solving the polynomial equation (32). Instead, to counteract roundoff-error problems, we solve the corresponding generalized eigenvalue problem as described in § 4.

### 3.5 Eigen-amplitude and maximum growth rate

The roots of (32) are denoted by , where . Since is singular when , can be determined by equation (30) up to a proportionality constant. Therefore, for each , the corresponding perturbed amplitude on the -th interface can also be determined up to a proportionality constant by equations (15) and (44). Let us call the vector , where we normalize , the “eigen-amplitudes” of frequency . Any arbitrary initial perturbation amplitudes can be decomposed into a linear combination of eigen-amplitudes, namely for some . Since the whole coupled system is linear, the small initial perturbation will grow as .

The growth rate for a single mode is since our time dependence is . If , then the mode is unstable. As described above, for an -layer system, we have different growth rates for a single , and we denote the largest growth rate by . Moreover, the maximum growth rate over all is denoted by , and the corresponding eigen-amplitudes of this most-unstable mode are denoted by ; denotes the corresponding wavenumber.

## 4 Generalized eigenvalue problem

It is well known that finding the roots of a polynomial via its coefficients is badly ill-conditioned Trefethen & Bau (1997). Correspondingly, we find that solving the determinant equation (32) directly, by treating it as a polynomial, is highly susceptible to roundoff errors when is not small. In particular, it is tempting to use the block structure of to reduce (32) to a determinant problem via a recurrence. However, the entries of this matrix are high-degree polynomials in whose coefficients thereby introduce roundoff ill-conditioning. Instead, we can turn this ill-conditioned root-finding problem into a generalized eigenvalue problem by exploiting the matrix-pencil structure of :

(33) |

Thus, finding the dispersion relation turns out to be a generalized eigenvalue problem with matrices . Since is non-singular, this (regular) generalized eigenvalue problem is typically well-conditioned Demmel & Kagstrom (1987) and can be solved via available numerical methods Anderson et al. (1999). In principle, further efficiency gains could be obtained by exploiting the sparsity of this pencil, but dense solvers are easily fast enough for up to hundreds.

## 5 Validation of our formulation

As a validation check, our -layer results can be checked against known analytical results in various special cases. We can also compare to previous finite-element calculations Deng et al. (2011).

### 5.1 Tomotika’s case:

Tomotika (1935) discussed the instability of one viscous cylindrical thread surrounded by another viscous fluid, which is equivalent to our model with . It is easy to verify that , where , gives the same determinant equation as (34) in Tomotika (1935).

In contrast with the Stokes-equations approach, Tomotika began with the full Navier–Stokes equations, treating the densities of the inner fluid and the outer fluid as small parameters and taking a limit to reach the Stokes regime. However, special procedures must be taken in order to obtain a meaningful determinant equation in this limit, because substituting and directly would result in dependent columns in the determinant. Tomotika proposed a procedure of expanding various functions in ascending powers of and , subtracting the leading terms in dependent columns, dividing a quantity proportional to , and finally taking a limit of and . We generalize this idea to the -shell problem in Appendix C, but such procedures are unnecessary if the Stokes equations are used from the beginning.

### 5.2 with equal viscosities

### 5.3 Navier–Stokes and inviscid cases

In Appendix C, we validate the generalized form of our instability analysis against previous results for inviscid and/or Navier–Stokes problems. For example, we find identical results to Chauhan et al. (2000) for the case of a two-fluid compound jet surrounded by air (with negligible air density and viscosity).

### 5.4 Comparison with numerical experiments

Deng et al. (2011) studied axisymmetric capillary instabilities of the concentric cylindrical shell problem for various viscosity contrasts by solving the full Navier–Stokes equation via finite-element methods. [The Stokes equation is a good approximation for their model, in which the Reynolds number is extremely low .] In particular, they input a fixed initial perturbation wavenumber , evolve the axisymmetric equations, and fit the short-time behaviour to an exponential in order to obtain a growth rate. With their parameters , , N/m, Pas, and , we compute the maximum growth rate for each ratio via the equation . For comparison, we also compute the growth rate for their fixed . [Because numerical noise and boundary artifacts in the simulations will excite unstable modes at , it is possible that and not will dominate in the simulations even at short times if the former is much larger.] The inset of figure 2 plots the wavenumber that results in the maximum growth rate versus the viscosity ratio (). In figure 2, we see that the growth rate obtained by Deng et al. (2011) (blue circles) agrees well with the growth rate predicted by linear stability analysis (blue line) except at large viscosity contrasts ( or ). These small discrepancies are due to the well-known numerical difficulties in accurately solving a problem with large discontinuities.

## 6 examples

In this section, we study the three-fluid () problem. Three or more concentric layers are increasingly common in novel fibre-drawing processes Hart et al. (2002); Kuriki et al. (2004); Pone et al. (2006); Abouraddy et al. (2007); Deng et al. (2008). By exploring a couple of interesting limiting cases, in terms of shell viscosity and shell thickness, we reveal strong connections between the case and the classic problem.

### 6.1 Case and : shell viscosity cladding viscosity

We first consider the limiting case in which the shell viscosity is much smaller than the cladding viscosities and . Substituting and into , equation (32) gives simple formulae for the growth rates

(34) |

and

(35) |

Note that is independent of , , and , while is independent of , , and . In particular, these growth rates are exactly the single-cylinder results predicted by Tomotika’s model, as if the inner and outer layers were entirely decoupled. This result is not entirely obvious, because even if the shell viscosity can be neglected, it is still incompressible and hence might be thought to couple the inner and outer interfaces. [Deng et al. (2011) conjectured a similar decoupling, but only in the form of a dimensional analysis.]

#### Case ,

In the regime that the shell viscosity is much smaller than the cladding viscosities and , we further consider the limit . It corresponds to the case with a high viscous fluid embedded in another low-viscosity fluid, which must of course correspond exactly to Tomotika’s case. From the asymptotic formulae of modified Bessel functions and for large arguments Abramowitz & Stegun (1992), we obtain

(36) |

Thus, the growth rate of possible unstable modes is given by in equation (34). Tomotika discussed this limiting case () and gave a formula (37) Tomotika (1935), which is exactly (34).

#### Case

The limit is equivalent to with a low-viscous fluid embedded in a high-viscosity fluid. For this case, it is easy to check that in equation (35) agrees with formula (36) in Tomotika (1935). However, we still have another unstable mode with a growth rate . Using the asymptotic formulae for and with small arguments Abramowitz & Stegun (1992), we obtain

(37) |

This extra unstable mode results from the instability of a viscous cylinder with infinitesimally small radius . In other words, (37) is the growth rate of a viscous cylinder in the air with a tiny but nonzero radius, which is also given by equation (35) of Lord Rayleigh (1892).

### 6.2 Thin shell case: ,

Next, we study a three-layer structure with a very thin middle shell; that is, with . A sketch of such a geometry is given in figure a. This is motivated by a number of experimental drawn-fibre devices, which use very thin (sub-micron) layers in shells hundreds of microns in diameter in order to exploit optical interference effects Hart et al. (2002); Kuriki et al. (2004); Pone et al. (2006).

Considering as a small parameter, we expand the determinant equation (32) in powers of . For a given wavenumber , the two roots of this equation are and , where can be computed analytically by dropping the terms of order in the determinant equation (32). After some algebraic manipulation, we find that actually is the solution for the structure (i.e., ignoring the thin shell) with a modified surface-tension coefficient . It is also interesting to consider a limit in which grows as shrinks. In this case, we find the same asymptotic results as long as grows more slowly than . Conversely, if it grows faster than , then the thin-shell fluid acts like a “hard wall” and all growth rates vanish. Instead of presenting a lengthy expression for , we demonstrate a numerical verification in figure 4. As indicated in figure a, the growth rate for approaches the the growth rate of with the summed surface-tension coefficients as . The parameters are and