2.1. Non-dimensionalisation and long-wave model

We proceed by non-dimensionalising the equations using the following scalings:

(2.18) \begin{equation} \begin{cases} x=Lx', \ r=r_{in}r',\ h_2=r_{in}h_2',\ h_1=r_{in}h_1', \ t=\frac {L}{u_{in}}t',\ {p=p_a+\frac {\mu _{hot}u_{in}}{L}p',} \\ u=u_{in}u',\ v=\frac {r_{in}u_{in}}{L}v',\ \mu (\theta )=\mu _{hot}\mu '(\theta ),\ \theta =\theta _{cold}+(\theta _{hot}-\theta _{cold})\theta ', \end{cases} \end{equation}

where primes denote the non-dimensional variables. After a straightforward substitution of these scalings into the governing equations and boundary conditions (2.1)–(2.17) we find that the system is described by the following dimensionless parameters:

(2.19) \begin{align} &\varepsilon =\dfrac {r_{in}}{L}, \qquad D= \frac {u_{out}}{u_{in}}, \qquad Re=\dfrac {\rho {u_{in}L}}{3\mu _{hot}}, \qquad P=\dfrac {p_HL}{\mu _{hot}{u_{in}}}, \qquad \Gamma =\dfrac {\gamma {L}}{\mu _{hot}{r_{in}u_{in}}}, \end{align}

(2.20) \begin{align} &Pe=\dfrac {\rho c_p u_{in}L}{k},\qquad b=B(\theta _{hot}-\theta _{cold}),\qquad G=\dfrac {\rho g L^2}{\mu _{hot}u_{in}},\qquad C=\dfrac {2h_wL}{\rho {c_p}r_{in}u_{in}}. \end{align}

Here, $\varepsilon$ is a measure of the slenderness of the fibre, $D$ is the draw ratio that represents the ratio of the output speed to the input speed, $Re$ is the Reynolds number which measures the ratio of inertial forces to viscous forces, $P$ is dimensionless internal hole pressurisation which measures the ratio of internal pressure forces in the hole to viscous forces, $\Gamma$ is the inverse capillary number which measures the ratio of surface tension to viscous forces, $Pe$ is the Péclet number which measures the ratio of heat advection to heat conduction, $b$ measures the relative change in the logarithm of the viscosity between the cold and hot temperatures, $G$ measures the ratio of gravitational to viscous forces and $C$ is a measure of the strength of the convective cooling relative to heat advection along the fibre.

The parameter values that occur in fibre drawing vary significantly between different applications. This is because drawing techniques are used for a broad range of fluids with very different material properties and for devices that have very different scales. A detailed description of parameter ranges that can occur can be found in Bechert (Reference Bechert2017). He showed that the largest value of $\varepsilon$ that occurs is $O(10^{-1})$ and it is often much smaller. In fibre drawing, $D$ is chosen to obtain the required area reduction from preform to fibre and is typically $O(10^3)-O(10^4)$ . However, there are other types of extrusional flows (Tronnolone et al. Reference Tronnolone, Stokes and Ebendorff-Heidepriem2017) for which such large values of $D$ are not appropriate and the effective draw ratio can be $O(1)$ . In applications, the value of $P$ will be zero if no internal pressure is applied. On the other hand, if internal pressure is applied, there will be a critical value of $P$ above which the hole will expand rapidly and burst (Chen et al. Reference Chen, Stokes, Buchak, Crowdy and Ebendorff-Heidepriem2015). The parameters $Re$ , $\Gamma$ , $Pe$ and $G$ can each take a broad range of both small and large values for various applications (Tronnolone et al. Reference Tronnolone, Stokes, Foo and Ebendorff-Heidepriem2016; Bechert Reference Bechert2017; Tronnolone et al. Reference Tronnolone, Stokes and Ebendorff-Heidepriem2017). In practice, the dimensionless cooling $C$ cannot be too large otherwise the viscosity of the fibre will become extremely large and enormous forces will be required to extend the fibre. The parameter $\alpha$ is the ratio of the hole radius to the fibre radius at the inlet. In applications this can vary widely since the optical properties of fibres differ over a very broad range of hole sizes. In extreme cases the area fraction of the glass in microstructured optical fibres may be less than 10 % (Argyros & Pla Reference Argyros and Pla2007). Due to this wide variation we perform a broad parameter study and focus on results that are unexpected or counter-intuitive.

Since the number of parameters is quite large and much of the behaviour is qualitatively similar to previous results for a solid fibre, we focus on the cases in which the presence of the hole in the fibre leads to surprising or counter-intuitive behaviour. In particular, we strongly focus on the role played by the relative size of the hole at the inlet, $\alpha$ , and the internal hole pressurisation, $P$ . With this in mind, we note that the role of axial thermal conduction in our system is to transport heat from the input aperture along the fibre. This effectively just reduces the effective cooling rate and thereby decreases the viscosity. Since this does not lead to any surprising behaviour, we consider the limit in which $Pe\gg 1$ and neglect the axial conduction terms. Since this represents a singular perturbation to the equations we also drop the boundary condition on the temperature at $x=1$ . In addition, the gravitational stretching term $G$ simply acts to modify the effective stretching rate along the axis of the fibre. This also does not lead to any new qualitative behaviour and hence we consider the limit $G\ll 1$ and neglect gravitation stretching. We note that these simplifications do not qualitatively change the interesting behaviour that we have found.

Following many previous studies (Yarin et al. Reference Yarin, Rusinov, Gospodinov and Radev1989; Dewynne et al. Reference Dewynne, Howell and Wilmott1994; Cummings & Howell Reference Cummings and Howell1999; Stokes et al. Reference Stokes, Tuck and Schwartz2000; Fitt et al. Reference Fitt, Furusawa, Monro, Please and Richardson2002; Wylie et al. Reference Wylie, Huang and Miura2007; He et al. Reference He, Wylie, Huang and Miura2016; Stokes et al. Reference Stokes, Wiley and Chen2019) we use the fact that the aspect ratio is a small parameter, that is $\varepsilon \ll 1$ , and expand all quantities in terms of $\varepsilon$ . The derivation of the equations satisfied by the leading-order quantities is straightforward, although relatively involved, and since it is contained in a number of previous papers we omit the details and refer the reader to Fitt et al. (Reference Fitt, Furusawa, Monro, Please and Richardson2002) and Stokes et al. (Reference Stokes, Wiley and Chen2019). Under the assumption that $\varepsilon ^2 Pe\ll 1$ and $\varepsilon ^2 Re\ll 1$ we obtain the following system of equations for the leading-order quantities:

(2.21) \begin{align} & Re\big(h_2^2-h_1^2\big)(u_t+uu_x)=[\mu (\theta )\big(h_2^2-h_1^2\big)u_x]_x+\Gamma [h_2+h_1]_x, \end{align}

(2.22) \begin{align} & \mu (\theta )\big(h_2^2-h_1^2\big)[(h_1^2)_t+(h_1^2u)_x]=P[h_1^2h_2^2]-\Gamma h_1h_2[h_2+h_1], \end{align}

(2.23) \begin{align} &\qquad\qquad\qquad \big(h_2^2-h_1^2\big)_t+[\big(h_2^2-h_1^2\big)u]_x=0, \end{align}

(2.24) \begin{align} &\qquad\qquad\qquad \big(h_2^2-h_1^2\big)\left (\theta _t+u\theta _x\right )= -C h_2\theta, \end{align}

where $\mu (\theta )=e^{-b\theta }$ . The leading-order boundary conditions are given by

(2.25) \begin{equation} u=1,\quad h_2=1,\quad h_1=\alpha, \quad \theta =1 \qquad \textrm {at} \quad x=0, \end{equation}

and

(2.26) \begin{equation} u=D \qquad \qquad \textrm{at} \qquad x=1. \end{equation}

Here, for convenience, we have dropped the primes that denote dimensionless variables. This system of equations has the same form as those derived by Fitt et al. (Reference Fitt, Furusawa, Monro, Please and Richardson2002) and will represent our main system of equations in this paper.

2.2. Steady-state solutions

To obtain the steady-state solutions, we set ${\partial }/{\partial t}\equiv 0$ in (2.21)–(2.24) and denote the dependent variables in the solutions, which no longer depend on time, with hats. Integrating the mass conservation equation (2.23) and using the boundary condition (2.25) at $x=0$ we obtain

(2.27) \begin{align} \hat {u}\big(\hat {h}_2^2-\hat {h}_1^2\big)&= (1-\alpha ^2). \end{align}

Using (2.27) to eliminate $\hat {u}$ from (2.21) and integrating with respect to $x$ we obtain

(2.28) \begin{align} Re\frac {(1-\alpha ^2)^2}{\hat {h}_2^2-\hat {h}_1^2} -\Gamma [\hat {h}_2+\hat {h}_1] +\frac {(1-\alpha ^2)\mu (\hat \theta )\big(\hat {h}_2^2-\hat {h}_1^2\big)_x}{\hat {h}_2^2-\hat {h}_1^2}={-F}, \end{align}

where $F$ is a constant that represents the pulling tension in the fibre. This constant must be chosen to satisfy the boundary conditions. Using (2.27) and (2.28) to simplify the remaining equations we obtain

(2.29) \begin{align} &\hat {h}_{1x}= \frac {-F\hat {h}_1}{2\mu (\hat \theta )(1-\alpha ^2)} -\frac {\Gamma \big(\hat {h}_2^2-\hat {h}_1^2\big)}{2\mu (\hat \theta )(1-\alpha ^2)}+\frac {P\hat {h}_1\hat {h}_2^2}{2\mu (\hat \theta )(1-\alpha ^2)} -\frac {Re(1-\alpha ^2)\hat {h}_1}{2\mu (\hat \theta )\big(\hat {h}_2^2-\hat {h}_1^2\big)}, \end{align}

(2.30) \begin{align} &\hat {h}_{2x}=\frac {-F\hat {h}_2}{2\mu (\hat {\theta })(1-\alpha ^2)}+\frac {\Gamma \big(\hat {h}_2^2-\hat {h}_1^2\big)}{2\mu (\hat \theta )(1-\alpha ^2)}+\frac {P\hat {h}_1^2\hat {h}_2}{2\mu (\hat \theta )(1-\alpha ^2)} -\frac {Re(1-\alpha ^2)\hat {h}_2}{2\mu (\hat \theta )\big(\hat {h}_2^2-\hat {h}_1^2\big)}, \end{align}

(2.31) \begin{align} &\qquad\qquad\qquad\qquad\qquad\qquad \hat \theta _x= -\frac {C}{(1-\alpha ^2)}\hat {h}_2\hat {\theta }. \end{align}

The boundary conditions at the inlet $x=0$ are given by

(2.32) \begin{equation} \hat {h}_1=\alpha, \quad \hat {h}_2=1, \quad \hat \theta =1. \end{equation}

Using (2.27) the boundary condition at the outlet $x=1$ is

(2.33) \begin{equation} \hat {h}_2^2-\hat {h}_1^2=(1-\alpha ^2)D^{-1} .\end{equation}

We note that the above system of equations is only valid if the hole size $\hat {h}_1\gt 0$ . If hole closure occurs ( $\hat {h}_1=0$ ) at a location $x=x_*$ in the domain, then for values of $x\gt x_*$ we must solve the system of equations appropriate for a solid fibre that are given by substituting $\hat {h}_1\equiv 0$ into (2.30) and (2.31).

In general, one must resort to numerical methods to solve the above system. However, there are some special cases in which we can obtain analytical results that shed significant light on some of the key mechanisms that govern the mechanics. In order to achieve this we make the following transformation:

(2.34) \begin{equation} A=\hat {h}_2^2-\hat {h}_1^2, \qquad H=\hat {h}_2+\hat {h}_1 \qquad \textrm {and} \quad x=(1-\alpha ^2)\int _0^y \mu (\hat {\theta }(\check {y})) {\rm d}\check {y}. \end{equation}

Here, $A$ is the dimensionless cross-sectional area, $H$ is the dimensionless total arclength in a cross-section and $y$ is a scaled axial distance. Under this transformation the system (2.29)–(2.33) simplifies dramatically and we obtain

(2.35) \begin{align} &\qquad\quad A_y= -F A + \Gamma H A -Re(1-\alpha ^2)^2, \end{align}

(2.36) \begin{align} & H_y=\frac {-FH}{2}+\frac {P(H^4-A^2)}{8H}-\frac {Re(1-\alpha ^2)^2 H}{2A}, \end{align}

(2.37) \begin{align} &\qquad\qquad\quad \hat \theta _y= -\frac {C(H^2+A)\hat {\theta }\mu (\hat {\theta })}{2H}. \end{align}

The boundary conditions at the inlet are given by

(2.38) \begin{equation} A=1-\alpha ^2,\quad H=1+\alpha, \quad \hat \theta =1 \qquad \textrm {at} \quad y=0. \end{equation}

Using (2.27) the boundary condition at the outlet is

(2.39) \begin{equation} A=(1-\alpha ^2)D^{-1}\qquad \textrm {at} \quad y=Y, \end{equation}

where $Y$ must be determined by solving $(1-\alpha ^2)\int _0^Y \mu (\hat {\theta }(\check {y}) ) {\rm d}\check {y}=1.$ These equations are only valid before hole closure occurs in which case $H^2\gt A$ . If hole closure occurs in the domain then we must solve (2.35) and (2.37) with $H^2\equiv A$ .

This transformation represents a simplification because $\hat \theta$ does not explicitly appear in (2.35) and (2.36). However, in order to determine $A$ and $H$ one needs to know the constant of integration $F$ and this must be found using the boundary condition (2.39) which does depend explicitly on $\hat \theta$ . Nevertheless, this transformation will help us to obtain some analytical expressions in some special cases that we will examine in detail.

2.3. Linear stability analysis

In order to determine the stability of the steady-state solutions there is no significant advantage in using the transformation to $A$ and $p$ and so we use the original variables. We add small perturbations of the form

(2.40) \begin{equation} \left. \begin{array}{r} u=\hat {u}(x)+\tilde {u}(x)e^{\lambda {t}},\\ h_1=\hat {h}_1(x)+\tilde {h}_1(x)e^{\lambda {t}},\\ h_2=\hat {h}_2(x)+\tilde {h}_2(x)e^{\lambda {t}},\\ \theta =\hat {\theta }(x)+\tilde {\theta }(x)e^{\lambda {t}}, \end{array}\right\} \end{equation}

where $\lambda$ is the growth rate of perturbations and tildes are used to represent the perturbation quantities. Substituting (2.40) into (2.21)–(2.24) and linearising in the perturbation quantities gives the following linear eigenvalue problem:

(2.41) \begin{align} & Re\big [\big(\hat {h}_2^2-\hat {h}_1^2\big)(\lambda {\tilde {u}}+\hat {u}\tilde {u}_x+\tilde {u}\hat {u}_x)+2(\hat {h}_2\tilde {h}_2-\hat {h}_1\tilde {h}_1)\hat {u}\hat {u}_x\big ]=\Gamma \big [\tilde {h}_2+\tilde {h}_1\big ]_x\nonumber \\ &\quad +\big [\mu (\hat \theta )\big(\hat {h}_2^2-\hat {h}_1^2\big)\tilde {u}_x+2\mu (\hat \theta )(\hat {h}_2\tilde {h}_2-\hat {h}_1\tilde {h}_1)\hat {u}_x-b\mu (\hat \theta )\tilde \theta \big(\hat {h}_2^2-\hat {h}_1^2\big)\hat {u}_x\big ]_x, \end{align}

(2.42) \begin{align} &\quad \mu {(\hat \theta )}\big(\hat {h}_2^2-\hat {h}_1^2\big)\big [2\lambda {\hat {h}_1\tilde {h}_1}+(\hat {h}_1^2\tilde {u}+2\hat {h}_1\tilde {h}_1\hat {u})_x\big ]\nonumber \\ &\quad +\big [2\mu {(\hat \theta )}(\hat {h}_2\tilde {h}_2-\hat {h}_1\tilde {h}_1) -b\mu {(\hat \theta )}\tilde {\theta }\big(\hat {h}_2^2-\hat {h}_1^2\big)\big ](\hat {h}_1^2\hat {u})_x \nonumber \\ &\quad =2P\big (\hat {h}_1^2\hat {h}_2\tilde {h}_2+\hat {h}_1\tilde {h}_1\hat {h}_2^2\big ) -\Gamma \big (\hat {h}_1^2\tilde {h}_2+2\hat {h}_1\tilde {h}_1\hat {h}_2+2\hat {h}_1\hat {h}_2\tilde {h}_2+\tilde {h}_1\hat {h}_2^2\big ), \end{align}

(2.43) \begin{align} &\qquad\quad 2\lambda \big (\hat {h}_2\tilde {h}_2-\hat {h}_1\tilde {h}_1\big )+\big [\big(\hat {h}_2^2-\hat {h}_1^2\big)\tilde {u}+2(\hat {h}_2\tilde {h}_2-\hat {h}_1\tilde {h}_1)\hat {u}\big ]_x=0, \end{align}

(2.44) \begin{align} &{\hat {u}\hat \theta _x\big (2\hat {h}_2\tilde {h}_2-2\hat {h}_1\tilde {h}_1\big )+\big(\hat {h}_2^2-\hat {h}_1^2\big) \big (\lambda \tilde \theta +\hat {u}\tilde \theta _x+\tilde {u}\hat \theta _x\big )=} -C\big [\hat {h}_2\tilde \theta +\tilde {h}_2\hat \theta \big ]. \end{align}

The boundary conditions are

(2.45) \begin{align} \tilde {u}=0,\quad \tilde {h}_1=0,\quad \tilde {h}_2=0, \quad \tilde {\theta }=0 &\qquad \textrm {at} \quad x=0, \nonumber \\ {\tilde {u}=0} &\qquad \textrm {at} \quad x=1. \end{align}

Equations (2.41)–(2.44) along with the boundary conditions (2.45) represent an eigenvalue problem for $\lambda$ . We note that these equations are only valid for steady-state solutions in which the hole does not close. However, in the study of stability, we will focus on cases in which hole closure does not occur.

3.1. Obtaining steady-state solutions

The system (2.29)–(2.31), or equivalently (2.35)–(2.37), represents three first-order ordinary differential equations with an undetermined constant of integration $F$ . There are three boundary conditions at the inlet and a condition at the outlet that allows us to determine the value of $F$ . As we will see below, in some special cases it is possible to obtain explicit expressions for the solutions using (2.35)–(2.37). However, in general, this is not possible and we must resort to numerical methods in which case we have solved (2.29)–(2.31). For a given value of $F$ it is straightforward to use the MATLAB routine ‘ode45’, which represents a fourth-order Runge–Kutta method, to accurately solve the system of ordinary differential equations (2.29)–(2.31). We employ this within a function whose input is $F$ and that outputs the cross-sectional area $\big(\hat {h}_2^2-\hat {h}_1^2\big)$ at $x=1$ . This forms the basis of a shooting technique in which one selects the value of $F$ to satisfy the boundary condition (2.33). The value of $F$ was determined using the MATLAB routine ‘fzero’ which uses a method based on bisection, secant and inverse quadratic interpolation techniques. The implementation of this method was very straightforward and obtaining the results did not pose any particular problems.

Figure 2. The effect of varying surface tension on the steady-state solution for a solid fibre and a hollow fibre is shown in columns (a) and (b), respectively. In column (a) the steady-state velocity $\hat {u}$ and outer radius $\hat {h}_2$ of a solid fibre are plotted as a function of $x$ for different values of the surface tension $\Gamma$ . In column (b) the same quantities are plotted along with the hole radius $\hat {h}_1$ for a fibre with $\alpha =1/2$ . In both cases the draw ratio $D=10$ , heat transfer coefficient $C=0$ , Reynolds number $Re=0$ and internal hole pressure $P=0$ .

3.2. Steady-state solutions for isothermal drawing

We begin by presenting some basic results that neglect inertia, pressure and cooling, and compare how surface tension affects the steady-state solutions for a solid fibre and a hollow fibre (figure 2). In the first column of figure 2 we see the results for the axial velocity $\hat {u}$ and outer radius $\hat {h}_2$ for a solid fibre for different values of the dimensionless surface tension $\Gamma$ . We note that this calculation is for a solid fibre and there is no profile for the radius of the hole $\hat {h}_1$ . The results for the different values of $\Gamma$ use the same draw ratio $D$ and so $\hat {u}$ and $\hat {h}_2$ adjust from $\hat {u}=\hat {h}_2=1$ at $x=0$ to $\hat {u}=D$ and $\hat {h}_2=D^{-1/2}$ at $x=1$ in all cases. The main effect of surface tension is to localise the thinning towards the outlet. The mechanism underlying this has been presented in Wylie et al. (Reference Wylie, Papri, Stokes and He2023). In fact, for a solid isothermal fibre with zero inertia an analytic solution may be obtained from (2.34), which becomes

(3.1) \begin{equation} A_y= -F A + \Gamma A^{3/2}. \end{equation}

This is a Bernoulli differential equation that can be solved subject to the boundary condition $A=1$ at $y=0$ to give

(3.2) \begin{equation} A=\left [\frac {\Gamma }{F}+\left (1- \frac {\Gamma }{F}\right )e^{Fy/2} \right ]^{-2}. \end{equation}

In the isothermal case with no hole, we have $\mu \equiv 1$ and $\alpha =0$ and so (2.39) becomes $A=D^{-1}$ at $y=1$ . Hence, we obtain an equation for $F$ given by

(3.3) \begin{equation} D^{1/2}=\frac {\Gamma }{F}+\left (1- \frac {\Gamma }{F}\right ) \exp \left (\frac {F}{2}\right ) . \end{equation}

For finite values of $\Gamma$ this equation can only be solved numerically. However, in the large $\Gamma$ limit, asymptotic calculations show that $F-\Gamma$ must be exponentially small. Hence, in this limit, (3.2) implies that there will be relatively little thinning over the bulk of the region, but there will be a layer of width $O(1/\Gamma )$ near the outlet where most of the thinning occurs.

In the second column of figure 2 we present the results for a hollow fibre with $\alpha =1/2$ . If the surface tension $\Gamma =0$ , then one can use (2.29) and (2.30) to readily see that $\hat {h}_2\hat {h}_{1x}-\hat {h}_1\hat {h}_{2x}=0$ . This implies that the ratio $\hat {h}_1/\hat {h}_2$ remains constant throughout the domain and so the hole will never close (Fitt et al. Reference Fitt, Furusawa, Monro, Please and Richardson2002). However, for sufficiently large surface tension the hole will close before the fibre reaches the outlet and $\hat {h}_1\equiv 0$ between the closure point and the outlet. The stronger the surface tension the more rapidly the hole closes. We note that for cases in which the hole closes, the conservation of mass dictates that the radius of the fibre at the outlet is given by $\hat {h}_2=(1-\alpha ^2)D^{-1}$ . On the other hand, if hole closure does not occur the outer radius of the fibre will be bigger than this value. In all cases, $\hat {u}$ adjusts from $\hat {u}=1$ at $x=0$ to $\hat {u}=D$ at $x=1$ and the increase in $\hat {u}$ becomes more localised towards $x=1$ as the surface tension increases. If the surface tension is sufficiently large that the hole closes then $\hat {h}_2$ adjusts from $\hat {h}_2=1$ at $x=0$ to $\hat {h}_2=(1-\alpha ^2)^{1/2}D^{-1/2}$ at $x=1$ and the decrease in $\hat {h}_2$ becomes more localised towards $x=1$ as the surface tension increases. However, if the surface tension is sufficiently small that the hole does not close, then increasing surface tension makes $\hat {h}_2$ uniformly smaller throughout the domain.

Figure 3. The effect of varying inertia on the steady-state solution for hollow fibres without surface tension ( $\Gamma =0$ ) and with surface tension ( $\Gamma =1$ ) is shown in columns (a) and (b), respectively. In both cases the draw ratio $D=10$ , inlet hole radius $\alpha =1/2$ , heat transfer coefficient $C=0$ , and internal hole pressure $P=0$ .

In figure 3 we examine how inertia affects the dynamics of hollow fibres in the absence of pressure and cooling; both with and without surface tension. In the first column of this figure we present the results for a hollow fibre with zero surface tension for different values of $Re$ and in the second column we present similar results when the surface tension is non-zero. As explained above, in this negligible surface tension case the ratio $\hat {h}_1/\hat {h}_2$ remains constant throughout the domain and so the hole can never close. Moreover, the variation in $\hat {u}$ , $\hat {h}_1$ and $\hat {h}_2$ becomes more localised towards $x=1$ as inertia is increased (see the first column of figure 3). This phenomenon was explained by Wylie et al. (Reference Wylie, Papri, Stokes and He2023). However, when the surface tension is non-zero the ratio $\hat {h}_1/\hat {h}_2$ is not constant and increasing $Re$ can cause the hole to close. This is shown in the second column of figure 3. In the presence of surface tension the variation in the velocity and the outer radius becomes more localised towards $x=1$ as the inertia increases. However, increased inertia makes the hole close more slowly with $x$ near the inlet and more rapidly with $x$ near the outlet.

Figure 4. The effect of varying the inlet hole radius on the steady-state solution for hollow fibres without surface tension ( $\Gamma =0$ ) and with surface tension ( $\Gamma =1$ ) is shown in columns (a) and (b), respectively. In both cases the draw ratio $D=10$ , heat transfer coefficient $C=0$ , Reynolds number $Re=0$ and internal hole pressure $P=0$ .

In figure 4 we examine how the size of the hole affects the dynamics of hollow fibres in the absence of inertia, pressure and cooling; both with and without surface tension. In the first column we present the results for a hollow fibre with zero surface tension for different values of $\alpha$ and in the second column we present similar results when the surface tension is non-zero. When the surface tension is zero the ratio $\hat {h}_1/\hat {h}_2$ remains constant throughout the domain and the velocity is independent of the hole size (see the first column of figure 4). In this case one can readily see that the outer radius $\hat {h}_2$ is also independent of the hole size if surface tension and pressure are zero. However, when the surface tension is non-zero, the dynamics is significantly more complicated. For small inlet hole radii the surface tension will cause the hole to close. For intermediate inlet hole radii the surface tension will not be strong enough to close the hole during the time that the fibre is in the domain. However, surprisingly, for inlet hole radii close to the inlet outer radius ( $\alpha$ close to unity) the hole will close very rapidly. There are some other interesting features of the dynamics. As $\alpha$ becomes close to unity and the hole closes very rapidly, the speed of the fibre decreases to conserve mass. There is also an associated decrease in the outer radius. After the hole closes the resulting fibre has to increase its speed from a low value up to $D$ . This is shown in the second column of figure 4. The rapid hole closure occurs because as $\alpha$ tends to unity the cross-sectional area becomes small and so offers little resistance to the surface tension force that remains finite. For negligible pressure, inertia and cooling, and assuming that hole closure does not occur, (2.35) and (2.36) can be solved explicitly to give further insight. After applying the boundary conditions (2.38) we obtain

(3.4) \begin{align} &\qquad\qquad\qquad H= (1+\alpha ) e^{-Fy/2}, \end{align}

(3.5) \begin{align} & A=(1-\alpha ^2)\exp \left [-Fy +\frac {2\Gamma (1+\alpha )}{F}(1-e^{-Fy/2}) \right ]. \end{align}

We note that this solution is only valid if $H^2\gt A$ , which corresponds to hole closure not occurring. In this case we can apply the boundary condition $A=(1-\alpha ^2)D^{-1}$ at $y=(1-\alpha ^2)^{-1}$ to obtain an equation for $F$ given by

(3.6) \begin{equation} \ln D=\frac {F}{1-\alpha ^2}-\frac {2\Gamma (1+\alpha )}{F}\left (1-\exp \left [-\frac {F}{2(1-\alpha ^2)}\right ]\right ). \end{equation}

This equation cannot be solved explicitly and we must use numerical methods to obtain $F$ . Having obtained $F$ , we can construct the solutions using (3.4) and (3.5).

If hole closure occurs at a location $y_*$ with $H=H_*$ and $A=A_*=H_*^2$ , then for values of $y\gt y_*$ we need to solve (3.1) subject to the boundary condition $A=H_*^2$ at $y=y_*$ . This can be solved explicitly to obtain

(3.7) \begin{equation} A=\left [\frac {\Gamma }{F}+\left (\frac {1}{p_*}- \frac {\Gamma }{F}\right )e^{F(y-y_*)/2} \right ]^{-2}. \end{equation}

Applying the boundary condition $A=(1-\alpha ^2)D^{-1}$ at $y=(1-\alpha ^2)^{-1}$ , we obtain an equation for $F$ given by

(3.8) \begin{equation} (1-\alpha ^2)^{-1/2}D^{1/2}=\frac {\Gamma }{F}+\left (\frac {1}{p_*}- \frac {\Gamma }{F}\right ) \exp \left (F\left [\left(1-\alpha ^2 \right)^{-1}-y_* \right ]\Bigg/2\right ) . \end{equation}

If $\alpha$ is close to one then one can readily show that $F-\Gamma H_*$ must be exponentially small. Using the approximation $F=\Gamma H_*$ , evaluating (3.4) and (3.5) at $y=y_*$ , and eliminating $y_*$ , we obtain an expression for $H_*$ given by

(3.9) \begin{equation} H_*=\frac {2(1+\alpha )}{2+\ln (1+\alpha )-\ln (1-\alpha )}. \end{equation}

We now reinterpret this solution in terms of the original variables for $\alpha$ close to one. The hole closure occurs on a scale $y\sim 1/F \sim 1/(H_*\Gamma )\,\sim -\ln (1-\alpha )/\Gamma$ as $\alpha \to 1$ . Since $x=(1-\alpha ^2)y$ , this corresponds to a scale of $x\,\sim -(1-\alpha )\ln (1-\alpha )/\Gamma$ which becomes small as $\alpha \to 1$ . This shows that there will be a narrow region near $x=0$ in which the hole rapidly closes. After the hole closes, we need to solve (3.1) with $F=\Gamma H_*$ and a boundary condition $A=H_*^2$ at $y=y_*$ . In the limit as $\alpha \to 1$ the surface tension forces become strong and the analysis is similar to that for a solid thread and large $\Gamma$ that we derived above. Therefore, after the hole closes most of the thinning will occur in a narrow region near the exit of width $y\sim -\ln (1-\alpha )/\Gamma$ or equivalently $x\,\sim -(1-\alpha )\ln (1-\alpha )/\Gamma$ . Hence, there will be very little thinning between the location of hole closure and a narrow region near the exit. This behaviour can clearly be observed in the curve for $\alpha =0.99$ in figure 4(b) that has an obvious plateau region after hole closure occurs. In the original variables, the height of the plateau can be obtained by recalling that $H=\hat {h}_2+\hat {h}_1$ and noting that at the point at which hole closure occurs we have $\hat {h}_2=H_*$ and $\hat {h}_1=0$ . Hence, the height of the plateau is simply $\hat {h}_2=H_*$ .

Figure 5. The hole size at the outlet, $\hat {h}_1(1)$ , is plotted against the hole size at the inlet, $\hat {h}_1(0)\equiv \alpha$ , for different values of surface tension $\Gamma$ for two different draw ratios, (a) $D = 5$ and (b) $D = 25$ . In both plots the heat transfer coefficient $C=0$ , internal hole pressure $P=0$ and Reynolds number $Re=0$ .

In figure 5 we plot $\hat {h}_1(1)$ , the size of the inner hole at the outlet, as a function of $\alpha$ , the size of the inner hole at the inlet, for different values of the surface tension $\Gamma$ . Since $A=(1-\alpha ^2)D^{-1}$ at the outlet, these results can be obtained using (3.6) and (3.4) evaluated at the outlet. The two plots are for two different values of the draw ratio $D$ . In both cases the inertia and hole pressurisation are negligible. For $\Gamma =P=0$ the ratio of $\hat {h}_1/\hat {h}_2$ is a constant throughout the entire length of the drawing. In this case, we can readily use (2.17) to see that $\hat {h}_1=\alpha D^{-1/2}$ at $x=1$ . That is, for $\Gamma =P=0$ the curve in figure 5 is a straight line. For small values of $\Gamma$ we see that both very small holes and very large holes will close, there is a range of hole sizes for which the hole does not close and there is a maximum achievable value of the hole size at the outlet. On the other hand, if the value of $\Gamma$ is sufficiently large then the hole will close no matter what the size of the hole at the inlet. This is consistent with the mechanism explained above since small inlet hole sizes close rapidly because the surface tension has very little work to do to close the hole, whereas large holes close rapidly because the surface tension induces very large hole closure rates when the hole radius approaches the outer radius. We note that the locations where the curves intersect with the horizontal axis correspond to hole closure occurring (i.e. $\hat {h}_1=0$ or $A=H^2$ ) exactly at $x=1$ or $y=(1-\alpha ^2)^{-1}$ . Evaluating (3.4) and (3.5) at $y=(1-\alpha ^2)^{-1}$ , using $A=H^2$ , and simplifying the result using (3.6) we obtain $F=(1-\alpha ^2)\ln [(1+\alpha )D/(1-\alpha ) ]$ . Substituting this back into (3.6) we obtain an expression for the critical value of $\Gamma$ above which hole close occurs

(3.10) \begin{equation} \Gamma _{crit}=\frac {(1-\alpha )(1+\alpha )^{1/2} \ln \left (\frac {1+\alpha }{1-\alpha }\right )\ln \left (\frac {(1+\alpha )D}{1-\alpha }\right )} {2\left [(1+\alpha )^{1/2}-(1-\alpha )^{1/2}D^{-1/2}\right ]}. \end{equation}

This expression is consistent with the numerical results shown in figure 5.

Figure 6. The steady-state solution for the velocity $\hat {u}$ , outer radius $\hat {h}_2$ and the inner radius $\hat {h}_1$ are plotted as a function of $x$ for different values of internal hole pressure $P$ and for surface tension $\Gamma =0$ in column (a) and $\Gamma =0.3$ in column (b). In both cases the draw ratio $D=10$ , heat transfer coefficient $C=0$ , Reynolds number $Re=0$ , and inlet hole radius $\alpha =1/2$ .

In figure 6 we examine how the size of the hole affects the steady states of hollow fibres in the absence of inertia and cooling, but include the effects of internal hole pressurisation and surface tension. In the first column of figure 6 we present the results for a hollow fibre with zero surface tension for different values of $P$ and in the second column we present similar results when the surface tension is non-zero. If $\Gamma =0$ we can use (2.28) to see that the cross-sectional area $\hat {h}_2^2-\hat {h}_1^2$ does not depend on $P$ . Using (2.27) we can, hence, deduce that $\hat {u}$ is also independent of $P$ if $\Gamma =0$ . This can be seen in the first column of figure 6. On the other hand, both $\hat {h}_1$ and $\hat {h}_2$ are affected by $P$ and both quantities uniformly increase as $P$ is increased. However, there is a finite value of $P$ at which both $\hat {h}_1\to \infty$ and $\hat {h}_2\to \infty$ in such a way that $\hat {h}_2^2-\hat {h}_1^2$ and hence $\hat {u}$ remains finite. This phenomenon was noted by Fitt et al. (Reference Fitt, Furusawa, Monro, Please and Richardson2002) and Chen et al. (Reference Chen, Stokes, Buchak, Crowdy and Ebendorff-Heidepriem2015), who described it as an ‘explosion’. For values of $P$ above the critical value no steady-state solutions exist.

In fact, the case of $Re=\Gamma =0$ is another situation in which the solution can be obtained analytically. Putting $\Gamma =Re=0$ in (2.35), integrating and applying the boundary condition (2.38) we obtain

(3.11) \begin{equation} A=\left(1-\alpha ^2\right)e^{-Fy}. \end{equation}

Substituting (3.11) into (2.36) and putting $Re=0$ we obtain

(3.12) \begin{equation} H_y=\frac {-FH}{2}+\frac {P\left (H^4-\left(1-\alpha ^2\right)^2e^{-2Fy}\right )}{8H}. \end{equation}

This equation can be made separable by considering the combination $H^2 e^{Fy}$ and after integrating and applying the boundary condition (2.38) we obtain

(3.13) \begin{equation} H=e^{-Fy/2}\left(1-\alpha ^2\right)^{1/2}\left (\frac {1+\alpha \exp \left [ \left(1-\alpha ^2\right)P\left(1-e^{-Fy}\right)/(2F)\right ]}{1-\alpha \exp \left [ \left(1-\alpha ^2\right)P\left(1-e^{-Fy}\right)/(2F)\right ]} \right )^{1/2}. \end{equation}

In the isothermal case, $\mu \equiv 1$ and applying (2.39) and (3.11) we obtain $F=(1-\alpha ^2)\ln D.$ An explosion will occur when the denominator in (3.13) becomes zero. The lowest value of $P$ for which this occurs can be obtained by setting $y=(1-\alpha ^2)^{-1}$ and yields

(3.14) \begin{equation} P_{crit}=-\frac {2 D \ln D \ln \alpha }{D-1}. \end{equation}

Using this expression for the parameters $\alpha =0.5$ , $D=10$ of figure 6 we see that the value at which an explosion first occurs for the case of zero surface tension ( $\Gamma =0$ , left column) is approximately 3.55. Results for $\Gamma \ne 0$ are shown in the second column of figure 6. In this case $\hat {u}$ depends weakly on $P$ . The quantities $\hat {h}_1$ and $\hat {h}_2$ also uniformly increase as $P$ is increased, but the surface tension acts to reduce the rate of increase with $P$ . The phenomenon of ‘explosion’ also occurs, but it happens at a larger value of $P$ than in the $\Gamma =0$ case.

3.3. Steady-state solutions for non-isothermal drawing

In figure 7 we examine how cooling affects the steady-state solution in the absence of inertia, surface tension and pressure. The first column shows the solution for different values of the heat transfer coefficient $C$ and a fixed value of the inlet hole radius $\alpha =1/2$ . As $C$ is increased the temperature decreases uniformly along the fibre. This temperature variation means that the viscosity near the inlet is lower than the viscosity near the outlet. Therefore, it is easier for the fibre to deform in the region near the inlet. This is reflected in the profiles for $\hat {u}$ , $\hat {h}_1$ and $\hat {h}_2$ . In the second column of figure 7 we examine how the inlet hole radius $\alpha$ affects the solution for a fixed value of the heat transfer coefficient $C=0.5$ . We immediately see that the temperature in the fibre decreases as the size of the hole increases. This can be easily understood by considering (2.24). The heat flux due to cooling to the environment is proportional to the outer radius, whereas the thermal mass (i.e. the amount of fluid to be heated) is proportional to $\hat {h}_2^2-\hat {h}_1^2$ . So, as the size of the inner hole, $\hat {h}_1$ , increases the convective flux is unaffected but the thermal mass decreases. This leads to the observed more rapid effective cooling as the size of the inner hole increases. The cooler temperatures for larger hole sizes mean that the viscosity increases more dramatically and the deformation becomes more localised near the inlet.

Figure 7. The steady-state solution for the temperature $\hat {\theta }$ , velocity $\hat {u}$ , outer radius $\hat {h}_2$ and inner radius $\hat {h}_1$ are plotted as a function of $x$ for (a) different values of heat transfer coefficient $C$ with inlet hole radius $\alpha =1/2$ and (b) different values of inlet hole radius $\alpha$ with heat transfer coefficient $C=0.5$ . In both cases the draw ratio $D=10$ , the viscosity parameter $b=6.3$ , Reynolds number $Re=0$ , surface tension $\Gamma =0$ and internal hole pressure $P=0$ .

Figure 8. The steady-state solution for the temperature $\hat {\theta }$ , velocity $\hat {u}$ , outer radius $\hat {h}_2$ and inner radius $\hat {h}_1$ are plotted as a function of $x$ for (a) different values of internal hole pressure $P$ with fixed heat transfer coefficient $C=0.5$ and inlet hole radius $\alpha =1/2$ , and (b) different values of heat transfer coefficient $C$ with fixed internal hole pressure $P=1$ and inlet hole radius $\alpha =0.95$ . In both cases the draw ratio $D=10$ , the viscosity parameter $b=6.3$ , Reynolds number $Re=0$ and surface tension $\Gamma =0$ .

We next examine how the interplay of cooling and internal hole pressure affects the steady-state profiles. The first column of figure 8 shows the solution for different values of the internal hole pressure $P$ with a fixed non-zero heat transfer coefficient and in the absence of inertia and surface tension. Here, we see that increasing $P$ acts to increase the effective cooling. This is because pressure (in the absence of surface tension) acts to increase the outer radius $\hat {h}_2$ and therefore the surface area over which heat is lost to the environment. As temperature decreases with $x$ the viscosity increases dramatically so that the deformation becomes more localised near the inlet, but the effect is relatively weak compared with the effect observed in figure 7. In the second column of figure 8 we examine how changing the cooling coefficient with a fixed non-zero pressure and in the absence of inertia and surface tension affects the steady-state solution. We have chosen a pressure and an inlet hole size sufficiently large that the hole would explode if the cooling rate were zero, so that we need to choose cooling coefficients in excess of a critical value in order to obtain solutions. The cooling acts to increase the viscosity, thereby preventing the pressure from causing an explosion of the inner hole. In the profiles for $\hat {h}_2$ we see that at small cooling coefficient $C\le 0.06$ the pressure has a dramatic effect in making the hole expand near the outlet. This larger outer radius means that the area over which cooling to the environment occurs is increased and so the temperature is smaller and the temperature gradient is larger near the outlet for smaller $C$ . Increasing $C$ above this range ( $C=0.07$ ) increases the cooling rate near the inlet and therefore the viscosity. Therefore the pressure has a weaker effect in expanding the hole leading to a smaller outer radius and smaller surface area through which heat is lost. This in turn results in weaker effective heat loss and an increase in the temperature at the outlet. However, as the cooling rate increases further ( $C\gt 0.07$ ) the temperature decreases rapidly and the viscosity increases dramatically. This makes deformation significantly more difficult, and hence pressure has a relatively small effect. Moreover, while the temperature gradient decreases towards the outlet, the temperature now decreases with increasing $C$ , due predominantly to the more rapid temperature drop near the inlet.

Figure 9. The temperature at the outlet $\hat {\theta }(1)$ plotted against the heat transfer coefficient $C$ for different values of inlet hole radius $\alpha$ . The draw ratio $D=10$ , internal hole pressure $P=1$ , surface tension $\Gamma =0$ , Reynolds number $Re=0$ and viscosity parameter $b=6.3$ .

In figure 9 we plot the temperature at the outlet as a function of the heat transfer coefficient for different initial hole sizes with pressure. For sufficiently small inlet hole sizes, the inner hole does not explode even in the case of no cooling ( $C=0$ ). For the parameters used in figure 9, this is true for values of $\alpha$ less than approximately 0.83. For these values of $\alpha$ the temperature at the outlet decreases monotonically with the heat transfer coefficient $C$ . However, if the inlet hole size is large enough and the heat transfer coefficient $C$ is low enough, the hole will explode and there will be no solution. For values of $C$ just above the critical value, the hole will expand very dramatically, and the heat flux will be large enough to almost cool the fibre down to the temperature of its surroundings by the time that it reaches the outlet. Further small increases in $C$ cause the hole to expand less dramatically and thus reduce the heat flux, so the temperature at the outlet increases. This continues until the temperature at the outlet reaches a maximum value, and for larger values of $C$ the temperature at the outlet decreases.

We note that the results presented in figures 7, 8 and 9 all have $\Gamma =Re=0$ . Therefore, the solutions (3.11) and (3.13) also hold for these cases. However, in this case, $C\ne 0$ , the viscosity $\mu$ will vary and so the $y$ location at which we have to apply the outlet boundary condition depends on the cooling rate. This means that one has to numerically solve a nonlinear equation to obtain the value of $F$ . This, in turn, implies that the profiles in the second column of figure 8 are not simple rescalings of one another as $C$ varies.

There are many other possible results for combinations of surface tension, inertia, cooling and hole pressure that we could present, for example, the interplay of the effects of surface tension and cooling. However, we have performed an extensive search of the parameter space and, in terms of steady states, all of the interesting and counter-intuitive behaviour that we have found has been detailed in the above figures.

4.1. Numerical methods for linear stability

For the most simple problem of a solid fibre with no hole, constant viscosity, no surface tension and no inertia, the most unstable eigenvalue is known to be complex (Matovich & Pearson, Reference Matovich and Pearson1969a ,Reference Matovich and Pearson b ). We therefore should also expect that, in general, the most unstable eigenvalue for our problem will be complex. A very accurate and efficient method to numerically determine eigenvalues of (2.41)–(2.44) with the boundary conditions (2.45) is to apply a shooting technique in which one supplements the boundary conditions $\tilde {u}=0$ , $\tilde {h}_1=0$ , $\tilde {h}_2=0$ and $\tilde {\theta }=0$ at $x=0$ with an additional condition $\tilde {u}_x=1$ at $x=0$ . The value of $\tilde {u}_x=1$ is arbitrary as a consequence of the linearity of the eigenvalue problem. Given a value of $\lambda$ , one can then easily solve the system (2.41)–(2.44) from $x=0$ to $x=1$ using the MATLAB routine ‘ode45’. For an arbitrary value of $\lambda$ , the solution of this computation will not satisfy the condition $\tilde {u}=0$ at $x=1$ . However, we can create a function whose input is $\lambda$ and that outputs the value of $\tilde {u}$ at $x=1$ . This will form the basis of a shooting technique in which one selects the value of $\lambda$ satisfying $\tilde {u}=0$ at $x=1$ . However, we note that in this case, the value of $\lambda$ will be a complex number. We therefore cannot use the MATLAB routine ‘fzero’ which only works for real-valued functions. Instead we use the MATLAB routine ‘fsolve’ which applies the Levenberg–Marquardt method and trust-region techniques. Although this method is fast and accurate it has the drawback that the operator may have an infinite number of eigenvalues and so it requires a good ‘initial guess’ for $\lambda$ to find the eigenvalue with the largest real part. One might naively imagine that one can overcome this difficulty by starting with a well-established value in a simple setting and using a continuation technique that slowly varies the values of the parameters until the desired result is obtained. However, it has been shown that for drawing problems with thermal effects the eigenvalue with the largest real part can change discontinuously with changes in the parameters (Wylie et al. Reference Wylie, Huang and Miura2007). We therefore have validated all of our results using a second-order finite difference method which converts (2.41)–(2.45) into a matrix eigenvalue problem. This method gives an approximation of the entire spectrum, but has the drawbacks that it is much slower and less accurate than the shooting method and also may give spurious results. Nevertheless, using the two techniques in combination we can confidently obtain the most unstable eigenvalue. In figure 10 we present a typical eigenvalue spectrum obtained using the finite difference method. The values of the most unstable eigenvalues can be accurately obtained using moderately small values of the grid spacing. However, the rapidly decaying modes are highly oscillatory and hence much smaller grid spacings are required to accurately resolve them. The spurious modes are always deep in the left half of the complex plane. Since we are only interested in the most unstable mode, it is straightforward to extract a good estimate of this mode from the finite difference method to use as an initial guess in the shooting method.

Figure 10. The most unstable few eigenvalues, plotted in the complex plane. These eigenvalues were obtained using a finite difference method. However, the stability results shown in the subsequent figures were obtained using a high accuracy shooting method that takes the most unstable eigenvalue given by the finite difference method as an ‘initial guess’. The parameters for generating this plot are $D=10$ , $\Gamma =1$ , $\alpha =1/2$ , $P=0$ , $Re=0$ and $C=0$ and the number of grid points was $N=400$ .

4.2. Stability of isothermal hollow fibres

We begin by considering the effect of the hole when surface tension is non-zero, but inertia, cooling and internal hole pressurisation are absent. Although this case was considered by Wylie et al. (Reference Wylie, Papri, Stokes and He2023), here, we present the results in a different format that will be useful for comparing with other results in which inertia, cooling and internal hole pressure are non-zero. Recall from sub§ 2.4 that the eigenvalue $\lambda$ indicates whether perturbations to a steady state grow ( $Re(\lambda )\gt 0$ ) or decay ( $Re(\lambda )\lt 0$ ) and hence shows the stability of the steady-state solution.

Figure 11. (a) The real part of the eigenvalue Re $(\lambda )$ versus the draw ratio $D$ for different values of the surface tension $\Gamma$ and inlet hole size $\alpha =1/2$ . (b) The critical draw ratio ( $D_{cr}$ ) plotted against the inlet hole size ( $\alpha$ ) for different values of surface tension $\Gamma$ . In both plots Reynolds number $Re=0$ , heat transfer coefficient $C=0$ and internal hole pressure $P=0$ .

Figure 11(a) shows the real part of the eigenvalue associated with the most unstable mode against the draw ratio $D$ for various values of the surface tension $\Gamma$ in the case that the inlet hole size is $\alpha = 1/2$ . For $\Gamma =0$ we see that draw ratios above a critical value of approximately $20.218$ are unstable, consistent with the classical result of Matovich & Pearson (Reference Matovich and Pearson1969b ). As $\Gamma$ increases the flow becomes less stable in the sense that instability occurs at lower values of the critical draw ratio. We also note that for sufficiently large values of $\Gamma$ the flow becomes unstable even for $D=1$ which indicates that no stable drawing to a slender fibre can occur. We denote the critical draw ratio by $D_{cr}$ , which has the property

(4.1) \begin{align} Re(\lambda ) \left \{ \begin{array}{ c c c} \lt 0& \textrm {for} & D\lt D_{cr},\\ =0& \textrm {for} & D=D_{cr},\\ \gt 0 & \textrm {for} & D\gt D_{cr}.\\ \end{array} \right . \end{align}

Figure 11(b) shows the critical value of the draw ratio $D_{cr}$ as a function of the inlet size of the hole $\alpha$ for various values $\Gamma \ge 0$ . We see that the role of the hole is always destabilising with $D_{cr}$ decreasing as $\alpha$ increases for any $\Gamma \gt 0$ .

Figure 12. The critical draw ratio ( $D_{cr}$ ) versus the inlet hole size ( $\alpha$ ) for different values of internal hole pressure $P$ for surface tension (a) $\Gamma =0$ and (b) $\Gamma =0.1$ . In both cases $C=Re=0$ .

Figure 13. The critical draw ratio ( $D_{cr}$ ) versus the inlet hole size ( $\alpha$ ) for different values of surface tension $\Gamma$ with $P=0.5$ , $Re=0$ and $C=0$ .

Figure 14. The critical draw ratio ( $D_{cr}$ ) versus the inlet hole size ( $\alpha$ ) for different values of Reynolds number $Re$ for surface tension (a) $\Gamma =0$ and (b) $\Gamma =0.1$ . In both cases $P=C=0$ .

Next we investigate the effect of adding internal hole pressure. In figure 12(a) we plot the results for $\Gamma =0$ . In the case of zero surface tension ( $\Gamma =0$ ) and no cooling ( $C=0$ ) the critical draw ratio remains independent of both the hole size $\alpha$ and the hole pressurisation $P$ . This is because if $\Gamma =0$ and $C=0$ then the temperature and hence the viscosity remain constant and (2.21) and (2.23) form a pair of equations for the quantities $\hat {u}$ and $\hat {h}_2^2-\hat {h}_1^2$ that decouple from (2.22). Since $P$ only appears in (2.22) and this is a first-order ordinary differential equation that advects information from the inlet towards the outlet, we can expect that the eigenvalue spectrum will be independent of $\alpha$ and $P$ if $\Gamma =0$ and $C=0$ . However, the structure is different for $\Gamma \gt 0$ . In figure 12(b) we plot $D_{cr}$ versus $\alpha$ for $\Gamma =0.1$ and different values of $P$ . In this case (2.21)–(2.23) are fully coupled and the stability depends on both $\alpha$ and $P$ . The main result is that both increasing the hole size or increasing the hole pressure act to destabilise the flow. This directly addresses the question raised by Fitt et al. (Reference Fitt, Furusawa, Monro, Please and Richardson2002). However, it is important to note that the destabilisation by pressure is only really significant for relatively large inlet hole sizes.

Figure 13 shows the critical draw ratio $D_{cr}$ as a function of hole size $\alpha$ for various values of the surface tension with hole pressurisation $P=0.5$ but in the absence of inertia and cooling. As explained above, the critical draw ratio is a constant in the case of no surface tension regardless of pressure. However, when the surface tension is non-zero, the critical draw ratio decreases with increasing inlet hole radius for $P\ge 0$ .

It is interesting to see if the stabilising effect of inertia can combine with the destabilising effect of surface tension to give a critical draw ratio that is a non-monotonic function of the inlet hole radius $\alpha$ . We consider this in figure 14, which gives results for $\Gamma =0$ and $\Gamma =0.1$ with $P=C=0$ . In the case of $\Gamma =0$ , the critical draw ratio is independent of $\alpha$ as seen in figure 14(a). As explained above, this is because (2.21) and (2.23) decouple from (2.22) and the equations of linear stability are identical to those for a solid fibre. This is also true in the case $P\ne 0, \Gamma =0, C=0$ . However, for $\Gamma \ne 0$ the results depend on $\alpha$ with an increase in $\alpha$ acting to destabilise the flow. Inertia still stabilises the flow, but its stabilising effect gets weaker as $\alpha$ increases. This shows that the stabilising effect of inertia cannot combine with the destabilising effect of surface tension to give non-monotonic behaviour with respect to $\alpha$ . In fact, this is true for other parameter combinations with inertia that we have examined. Therefore, in what follows, we will only consider $Re=0$ and will focus on combinations of parameters that give more complicated and unexpected non-monotonic behaviour.

Figure 15. The critical draw ratio ( $D_{cr}$ ) versus the inlet hole size ( $\alpha$ ) for different values of heat transfer coefficient $C$ with $P=0, Re=0,\Gamma =0$ and $b=6.3$ .

Figure 16. The critical draw ratio ( $D_{cr}$ ) versus the inlet hole size ( $\alpha$ ) (a) for different values of surface tension $\Gamma$ with heat transfer coefficient $C=0.1$ , (b) for different values of heat transfer coefficient $C$ with $\Gamma =0.1$ . In both cases $P=0$ , $Re=0$ and $b=6.3$ .

4.3. Stability of non-isothermal hollow fibres

In figure 15 we show how the heat transfer coefficient $C$ affects stability when $Re=\Gamma =P=0$ . As already seen, if $C=0$ the critical draw ratio is independent of the inlet hole size $\alpha$ . However, the heat flux term in (2.24) is proportional to ${h}_2$ and the viscosity in (2.21) depends on $\theta$ , so if $C\ne 0$ then all four equations (2.21)–(2.24) are fully coupled. We see that increasing $C$ acts to increase $D_{cr}$ and hence stabilise the flow. This is the same effect as found for the case of a solid fibre by Shah & Pearson (Reference Shah and Pearson1972a ). Furthermore, increasing $\alpha$ also acts to stabilise the flow. This can be understood in terms of the mechanism that we explained in the steady-state results, namely an increase in the size of the hole decreases the thermal mass while leaving the heat flux term unchanged which increases the rate at which the temperature decreases, so stabilising the drawing of hollow fibres.

Figure 16 shows how the interplay of cooling and surface tension affects stability for drawing of hollow fibres with different inlet hole sizes in the absence of pressure and inertia. In figure 16(a), we see that in the case of zero surface tension the critical draw ratio increases monotonically with increasing $\alpha$ . This stabilisation is as a result of the enhanced cooling effect described above. On the other hand, for non-zero surface tension, we see that increasing $\alpha$ acts to destabilise for most values of $\alpha$ , but acts to stabilise at very large values of $\alpha$ . This is because, for small values of $\alpha$ , the effect of increasing $\alpha$ on the thermal mass of the fibre is rather weak (the variation is quadratic in $\alpha$ ) and thus the enhanced cooling and consequent thermal stabilisation are rather weak. On the other hand, increasing $\alpha$ results in a more pronounced increase (linear in $\alpha$ ) in the surface area over which surface tension acts. Therefore, the destabilising effect of surface tension dominates. For values of $\alpha$ close to unity the surface area approaches a constant whereas the thermal mass becomes very small. This means that the thermal stabilisation effect will dominate the destabilising effect of surface tension. This interplay is also seen in figure 16(b). The critical draw ratio monotonically decreases with increasing inlet hole size in the absence of cooling because of the destabilising effect of surface tension. If the cooling is sufficiently small, this situation remains, excepting for very large $\alpha$ when the small thermal mass and fast cooling stabilises the flow. However, if cooling is sufficiently large, then its stabilising effect can dominate for all values of $\alpha$ .

Figure 17. The critical draw ratio ( $D_{cr}$ ) versus the inlet hole size ( $\alpha$ ) for different values of internal hole pressure $P$ for heat transfer coefficient $C=0.05$ , $\Gamma =0$ , $Re=0$ and $b=6.3$ . Panel (a) shows a close-up of the small- $\alpha$ behaviour of (b).

Figure 18. The critical draw ratio ( $D_{cr}$ ) versus the inlet hole size ( $\alpha$ ) for different values of heat transfer coefficient $C$ with $P=1, Re=0,\Gamma =0$ and $b=6.3$ .

In figure 17 we show how the interplay between internal hole pressurisation and cooling affects the stability for different inlet hole sizes when inertia and surface tension are zero. We recall that for a given pressure there is a maximum value of $\alpha$ above which the hole will explode. This is the reason that the curves terminate at values $\alpha \lt 1$ . Figure 17(a,b) gives results for weak cooling. Recall that in figure 12(a) we saw that the critical draw ratio is not affected by internal pressurisation or the inlet hole size if the cooling coefficient $C$ and surface tension $\Gamma$ are zero. Similarly, for weak cooling the graphs show that for small inlet hole sizes the presence of the hole has a very weak effect and the critical draw ratio is almost independent of $\alpha$ . In figure 17(a,b) we see that for weak cooling the flow is stabilised with increasing $\alpha$ up to a critical value of $\alpha$ which depends on the pressure, beyond which further increases in $\alpha$ act to destabilise. This is because cooling is the dominant effect for small values of $\alpha$ and this acts to stabilise the flow while for large values of $\alpha$ the destabilising effect of internal hole pressure is dominant. We also see that internal hole pressurisation can play both a stabilising and destabilising role. On the one hand, for values of $\alpha$ below the critical value the internal hole pressurisation stabilises as seen in figure 17(b). This is because it acts to increase the outer radius of the fibre, enhances the effective cooling rate and hence stabilises. On the other hand, for values of $\alpha$ above the critical value the internal hole pressurisation destabilises as a result of the mechanism that we observed in the isothermal case.

In figure 18 we further examine how the interplay between cooling and hole pressurisation affects the stability for different inlet hole sizes when inertia and surface tension are zero. We immediately see that if $C=0$ the critical draw ratio is independent of the hole size and is given by the classical value of approximately 20.218 that corresponds to a solid fibre with no inertia and surface tension. This was explained above in discussing figure 15. For small values of $\alpha$ the critical draw ratio increases with increasing $\alpha$ because of reduced thermal mass and faster cooling as seen for $P=0$ in figure 15. However, as $\alpha$ approaches the value at which explosion occurs there is a competition between the stabilising effect of cooling and the destabilising effect of internal pressure. For large values of $\alpha$ and small to moderate values of $C$ the destabilising effect of internal pressure dominates and the critical draw ratio decreases rapidly with increasing $\alpha$ . For larger values of $C$ the stabilising effect of cooling wins and the critical draw ratio is an increasing function of $\alpha$ over the whole range of $\alpha$ .

Figure 19. The critical draw ratio ( $D_{cr}$ ) versus the inlet hole size ( $\alpha$ ) (a) for different values of heat transfer coefficient $C$ with $\Gamma =0.03$ , (b) for different values of surface tension $\Gamma$ with $C=0.05$ . In both cases $P=1$ , $Re=0$ and $b=6.3$ .

Figure 19(a) shows how small and non-zero surface tension ( $\Gamma =0.03$ ) changes the results shown in figure 18 for which surface tension was zero. In the zero surface tension case, the critical draw ratio was independent of $\alpha$ for $C=0$ , but as seen in figure 19(a), the presence of surface tension destabilises the flow and the critical draw ratio is a monotonic decreasing function of $\alpha$ . When cooling is non-zero, surface tension continues to destabilise but cooling stabilises, therefore figure 19(a) shows similar qualitative behaviour to that shown in figure 18. In figure 19(b) we fix the cooling coefficient ( $C=0.05$ ) and show how surface tension affects the stability. For very small values of $\alpha$ surface tension is the dominant effect and this acts to destabilise the flow. This effect becomes weaker as $C$ increases or $\Gamma$ decreases, but for all values of $\Gamma \ne 0$ the destabilising effect can still be observed even though it is very localised near $\alpha =0$ for large $C$ or small $\Gamma$ . For small surface tension the interplay between the stabilising effects of cooling and the destabilising effects of internal pressure mean that the critical draw ratio increases and then decreases as $\alpha$ is increased. However, if the surface tension is large enough its destabilising effect will overcome the stabilising effect of cooling and so the critical draw ratio is a monotonic decreasing function of $\alpha$ .

In summary, we have shown that there is a complicated interplay between the destabilising effects of internal pressure and surface tension and the stabilising effects of cooling and inertia. As the inlet hole size $\alpha$ increases the stabilising effect of cooling increases and the destabilising effects of surface tension and pressure increase. Internal pressure has a weak effect excepting near $\alpha =1$ when the fibre is close to ‘exploding’(or bursting). It is ‘then’ the interplay between surface tension and cooling that are most important. The destabilising effect of surface tension increases with $\alpha$ because the amount of surface area increases. At the inlet the combined arclengths of the inner and outer surfaces are $2\pi (1+\alpha )$ and this gives an indication of how important the effect of surface tension can be near the inlet. In the limit $\alpha \to 1$ the arclength near the inlet is double the value for $\alpha =0$ . The stabilising effect of cooling has a very different relationship with $\alpha$ . As the hole size increases the thermal mass of the hollow fibre is proportional to $\hat {h}_2^2-\hat {h}_1^2$ . This quadratic dependence on the inner hole radius means that small holes have a very small thermal stabilising effect. However, as $\alpha \to 1$ the thermal mass tends to zero and so the temperature of the fibre will very rapidly decrease to the ambient temperature. This means that the thermal stabilising effect will be very strong for large values of $\alpha$ . Nevertheless, since the temperature must ultimately saturate we expect that there will be a limit to the largest thermal stabilising effect that we can achieve. The weak destabilising effect of internal pressure for small $\alpha$ is also due to the small (internal) boundary length on which it acts. It is difficult to assess how large the destabilising effects of internal pressure are in the limit as the explosion is approached. Nevertheless, our results indicate that the stabilising effects of cooling can overcome the destabilising effects of internal pressure even very close to explosion. Since, it is clear that the thermal stabilising effect is bounded, this indicates that the destabilising effects of internal pressure must also be bounded.

Figure 20. The critical draw ratio ( $D_{cr}$ ) versus the inlet hole size ( $\alpha$ ) for large values of heat transfer coefficient $C$ with large $P=2$ , large $\Gamma =1$ , $Re=0$ and $b=6.3$ .

Figure 20 illustrates how sufficiently large cooling coefficient influences the critical draw ratio for large pressure and large surface tension in the absence of inertia. This graph clearly shows there are two distinct effects. Because of the interaction of sufficiently large surface tension with pressure, the critical draw ratio decreases as inlet hole size $\alpha$ increases, for small to moderate $\alpha$ . However, when $\alpha$ becomes sufficiently large, the action of sufficiently large cooling coefficient due to large $C$ and expansion of the hole under large pressure resulting in increased external boundary length, brings a rapid increase in the critical draw ratio as $\alpha$ increases.