Asymptotic Modeling of Optical Fibres: Annular Capillaries and MOFs

1. Introduction

Microstructured optical fibres (MOFs) also known as Photonic Crystal Fibres (PCFs) or holey fibres are a novel kind of optical fibres that arose approximately thirty years ago and led to a revolution in the field of optical fibres technology. These fibres are notable for their ability to manipulate, guide, and control light, enabling a wide range of applications in different fields [1,2]. They contain a highly controlled and periodic array of air holes or voids of micron-scale size extending through the whole fibre length that is precisely designed and structured within the fibre. Light guidance within MOFs mainly occurs by two mechanisms, i.e., the index guiding or the photonic bandgap (PBG) effect. The former is achieved by replacing one or more air capillaries in the center of the whole structure with solid glass rods, and it is essentially a modified form of the total internal reflection (TIR) mechanism since the resulting solid core region possesses a higher refractive index than the surrounding region laced with air holes. On the other hand, by removing one or more air capillaries from the centre of the air hole array or by replacing them with one air capillary of different size, the PBG effect is obtained [3,4]. MOFs are highly versatile and flexible and permit to obtain a wide range of optical properties that are not possible to be realized with conventional optical fibres. This is possible since the arrangement and the sizes of the air holes and the crystal lattice constant can be suitably adjusted to select the guidance of specific modes. PCFs are widely employed in high-power and energy transmission [5,6], fibre amplifiers [7], Brillouin scattering [8], supercontinuum generation [9], optical communications [10], stimulated Raman scattering [11], and optical sensors [12,13], among others.

The first step of the PCFs manufacturing process is to create definite structures called preforms, which are macroscopic versions of the final fibre structures containing periodic patterns of air holes in their cross-section. A popular and versatile technique utilised to create preforms is the stack-and-draw method, where drawn capillary tubes and solid rods are precisely aligned, stacked, and fastened together. This method allows structures with a complex pattern of air holes to be easily, flexibly and quickly manufactured. Subsequently, preforms are heated and drawn in high-tech furnaces into long and thin fibres known as canes or preform-canes. These structures do not attain the required final fibre dimensions and one or more drawing stages are required to obtain the final size. The hole dimensions and positions within a fibre structure are of paramount importance to obtain a device with specific optical properties, and only small variations concerning the hole size and locations are allowed. Deviations from the designed structure may lead to a fibre with undesired optical characteristics. The alteration of the shape and the size of holes during the drawing process or even their collapse may occur due to the complex interplay between the surface tension effects and the stress redistribution within the fibre. Often, internal pressurization is utilized to avoid hole collapse and to precisely control the size of the holes. In general, the fabrication of PCFs is a complex process that requires great ability, precision, control, and expertise to preserve the required highly regular fibre structure from the beginning to the end of the drawing. Experimentally, only a few parameters can be controlled, that is, the peak temperature of the furnace, the feed and drawing speed, and the value of the internal pressurization. Usually, many preform-canes are utilized for the development of new fibre structures and for process optimization. However, this experimental trial-and-error routine is costly and time-consuming. This burden motivated the development of many mathematical models and utilization of numerical simulations that describe the drawing process over the years since they can potentially reduce the experimental costs and improve the physical understanding of the whole process.

Modelling of drawing of axis-symmetric molten treads was first driven by the spinning of synthetic fibres for textiles. To this end, Matovich and Pearson [14], Peak and Runk [15], Glicksman [16], and Myers [17] utilized a perturbation approach to derive one-dimensional models that included asymptotic mass, momentum, and energy equations with an increasing degree of sophistication. Modelling of drawing of axis-symmetric optical fibres have been undertaken by Yarin [18], Fitt et al. [19], and Luzi [20,21], among others. They derived one-dimensional models using an asymptotic expansion in terms of the small fibre aspect ratio. They investigated the effects of inertia, gravity, surface tension, and inner pressurization in the process of fibre drawing and demonstrated that an appropriate choice of furnace temperature distribution permits to obtain a good concordance between experimental trials and numerical results. Subsequently, Voyce et al. [22,23] integrated the effect of fibre rotation into the model of Fitt et al. [19]. Preform rotation is a very useful means to control fibre size and optical effects. More recently, Taroni et al. [24], and Luzi et al. [25] extended the work by Fitt et al. [19] by deriving asymptotic energy equations both for a solid fibre and a capillary to investigate the relevance of conduction, convection, and radiation phenomena in the drawing process. All these models mentioned so far have been derived in an Eulerian reference frame, where fibre drawing is a steady-state problem. Despite the fact that the formulation of the problem in the Eulerian coordinate system greatly simplifies the original full three-dimensional problem in the case of axis-symmetric geometries, it is not suitable for non-axis-symmetric geometries.

On the other hand, the utilization of a Lagrangian coordinate system with an asymptotic approach allows for modeling fibres with general cross-sectional shapes. Works in this direction have been initiated by Dewinne et al. [26,27] who modelled the drawing of a solid viscous thread with an arbitrary cross-sectional shape. In their works, they neglected the presence of surface tension but they examined the effects of inertia and gravity. Griffiths and Howell [28,29] investigated the evolution of a thin-annular two-dimensional fluid region retaining the effects of surface tension and applied pressure. Cummings and Howell [30], Stokes et al. [31], and Chen et al. [32] developed an asymptotic formulation for modelling the fibre drawing process. To this end, they employed a Lagrangian reference frame with asymptotic techniques and a reduced time variable $\tau$. By doing so, they showed that the full three-dimensional problem can be split into a one-dimensional axial-stretching problem and a two-dimensional Stokes flow free-boundary problem driven by unit surface tension and unit viscosity. Both problems are coupled by the whole length of the cross-sectional perimeter. This formulation is very general and applies to fibres with any cross-sectional shape and variable viscosity. Besides, it allows for including the effects of the inner pressurization as well. Subsequently, Stokes et al. [33] derived a three-dimensional asymptotic energy equation and included it in their previous models. To simplify the treatment of the cross-sectional problem, Buchak et al. [34] proposed the Generalized Elliptical Pore Model (GEPM) to model the drawing process of fibres containing many holes. This model is based on the complex variable formulation of the cross-plane problem and the utilization of inner-outer-matched asymptotic expansions. Assuming an elliptical shape of the holes that compose the fibre cross-section, their evolution is dictated by the solution of a set of nonlinear ordinary integrodifferential and algebraic equations. The GEPM is an extension of the model developed by Crowdy [35] to investigate the evolution of a two-dimensional inviscid elliptical bubble. In the case of fibres with arbitrary cross-sectional hole shapes, the numerical scheme proposed by Buchak and Crowdy [36] allows for a precise representation of evolving free surfaces of arbitrary shape even with extreme curvature gradients.

Alternatively, the drawing process of optical fibres with an arbitrary cross-sectional shape can also be simulated using the finite-element method (FEM). To this end, Xue et al. [37,38,39], and Luzi et al. [40] carried out numerical computations of the drawing process of optical fibres with a handful of holes, examining how the hole shape changes throughout the drawing due to the complex interplay among surface tension, inner pressure, and redistribution of viscous stresses. Besides, Luzi et al. [40] showed that a good agreement between experiments and numerical simulations can be achieved by specifying a suitable furnace temperature profile. However, three-dimensional FEM simulations are computationally expensive, and high-performance computers (HPCs) are usually required. On the contrary, Chen et al. [41] showed that the coupling between asymptotic techniques and the numerical method proposed by Buchak and Crowdy [36] leads to an enormous reduction of the computational time compared to the FEM simulations. Besides, the agreements between numerical results and experimental trials are superior to those obtained using the FEM. This demonstrates that the asymptotic approach constitutes a powerful modelling tool for describing the fibre drawing process.

In this manuscript, we model the drawing process of annular capillaries and holey fibres both in the absence and in the presence of internal pressurization utilising the asymptotic models of Stokes et al. [31], Chen et al. [41], and Buchak et al. [34]. To this end, we do not model the thermal effects due to the heat transfer during the drawing process, but we include them in the viscosity model by specifying a suitable functional form of the fibre temperature profile. The assumption of thermal equilibrium between the furnace walls and the fibre during the drawing process is only realistic for fibre velocities of the order of magnitudes of the feeding speeds [23,25]. Nevertheless, we demonstrate here that a good accordance between experimental and numerical results can be achieved by a suitable specification of the fibre temperature profile. The remainder of the manuscript is organized as follows. In Section 2, we expose the theoretical formulation of the problem, providing a concise description of the mass, momentum, and of corresponding boundary conditions that govern the drawing process of annular capillaries and holey fibres. In Section 3, we compare the numerical results with experimental data available in the literature, and in Section 4, we discuss the results, the limitations of the present approach, and hypothesize modalities for its improvement.

2. Model Description

2.1. Three-dimensional model

We consider the problem of a fibre preform that enters the top inlet of a furnace with a constant feeding speed $U_0$. The fibre is heated inside the oven and is pulled from the bottom so that it reaches a final fixed speed $U_F$ at the exit of the furnace. The fibre cross-section has many air holes which greatly shrink during the drawing process and gradually attain the required final size at the outlet of the oven. The drawing process is depicted in Figure 1. To model the drawing process of PCFs, we start from the three-dimensional mass and momentum equations written in the Cartesian coordinate system and we assume an incompressible Newtonian fluid [42]

$$\nabla ·\mathit{u}=0,$$

(1)

$$\begin{array}{c}\hfill \rho \left(\frac{\partial \mathit{u}}{\partial t}+\mathit{u}·\nabla \mathit{u}\right)=-\nabla p+\nabla ·\sigma \end{array}$$

(2)

The origin $x=y=z=0$ of the coordinate system is located in the middle of the inlet of the furnace and in the geometrical centre of the fibre cross-sectional area at the entrance of the oven. y and z are the coordinates in the cross-sectional plane and x is the coordinate in the drawing direction. $\mathit{u}=(u,v,w)$ and p denote the velocity vector and the pressure of the fluid. t, $\rho$ and $\sigma$ are the time, density, and the stress tensor of the fluid, respectively. $\sigma$ reads

$$\sigma =\mu \left(T\right)\left(\nabla \mathit{u}+{\left(\nabla \mathit{u}\right)}^T\right),$$

(3)

where $\mu \left(T\right)$ represents the dynamic viscosity of the fluid. The two terms on the left-hand side of Eq. (2) denote the temporal and convective inertial acceleration. The three terms on the right-hand side represent the pressure gradient, the divergence of the stress tensor, and the sum of the body forces. The kinematic and stress conditions read

$$\begin{array}{cccc}\hfill \frac{\partial G^{\left(i\right)}}{\partial t}+\mathit{u}·\nabla G^{\left(i\right)}& =0,\hfill & \hfill i& =0,1,....,N,\hfill \end{array}$$

(4a)

$$\begin{array}{cccc}\hfill \sigma ·{\mathit{n}}^{\left(i\right)}& =-\left(\gamma {\kappa }^{\left(i\right)}+p_H^{\left(i\right)}\right){\mathit{n}}^{\left(i\right)},\hfill & \hfill i& =0,1,....,N,\hfill \end{array}$$

(4b)

where $G^{\left(0\right)}(x,y,z,t)=0$ denote the external free surface of the fibre and $G^{\left(i\right)}(x,y,z,t)=0$, $i=1,2,...,N$ indicate the free surface of each of the internal channel of the fibre. $p_H^{\left(i\right)}$ and ${\kappa }^{\left(i\right)}$ are the applied internal pressurization and the local curvature of the ith free surface. Besides, the operator ∇ in Cartesian coordinates reads

$$\nabla =\frac{\partial }{\partial x}\widehat{\mathit{i}}+\frac{\partial }{\partial y}\widehat{\mathit{j}}+\frac{\partial }{\partial z}\widehat{\mathit{k}}$$

where $\widehat{\mathit{i}}$, $\widehat{\mathit{j}}$, and $\widehat{\mathit{k}}$ denote the unit vectors in the x, y, and z direction, respectively. The unit normal vector ${\mathit{n}}^{\left(i\right)}$, $i=0,1,....,N$, may be written as

$${\mathit{n}}^{\left(i\right)}=\frac{\nabla G^{\left(i\right)}}{|\nabla G^{\left(i\right)}|}$$

2.2. Non-Dimensionalization

Cumming and Howell [30], Stokes et al.[31], and Chen et al. [32] suggest the following scalings

$$\begin{array}{cccccc}\hfill (u,v,w)& =U_0(\overline{u},ϵ\overline{v},ϵ\overline{w}),\hfill & \hfill t& =\frac{L}{U_0}\overline{t},\hfill & \hfill p& =\frac{{\mu }_s{U}_0}{L}\overline{p},\hfill \\ \hfill (x,y,z)& =L(\overline{x},ϵ\overline{y},ϵ\overline{z}),\hfill & \hfill \mu \left(x\right)& ={\mu }_s\overline{\mu }\left(\overline{x}\right),\hfill & \hfill {\kappa }^{\left(i\right)}& =\frac{1}{ϵL}{\overline{\kappa }}^{\left(i\right)},\hfill \\ \hfill \chi & ={\chi }_0\overline{\chi },\hfill & \hfill \mathrm{\Gamma }& =ϵL\overline{\mathrm{\Gamma }},\hfill & \hfill \gamma & =\frac{{\mu }_s{U}_0{\chi }_0}{L}\overline{\gamma },\hfill \end{array}$$

to non-dimensionalize the mass and momentum equations and the kinematic and dynamic boundary conditions. L represents the length of a “hot zone" in the furnace where the viscosity of the glass is low enough to enable the drawing. ${\mu }_s$ denotes the glass viscosity at the softening point, ${\chi }_0$ and $\chi$ are the square of the cross-sectional area of the fibre at the entrance of the furnace and at any location in the drawing direction, respectively. $\mathrm{\Gamma }$ represents the total sum of the perimeters of the outer and inner cross-sectional free surfaces. Chen et al. [32] show that using this scaling the mass and momentum equations assume the form

$$\frac{\partial \overline{u}}{\partial \overline{x}}+\frac{\partial \overline{v}}{\partial \overline{y}}+\frac{\partial \overline{w}}{\partial \overline{z}}=0$$

(5)

$$\begin{array}{c}\hfill Re{ϵ}^2\left[\frac{\partial \overline{u}}{\partial \overline{t}}+\overline{u}\frac{\partial \overline{u}}{\partial \overline{x}}+\overline{v}\frac{\partial \overline{u}}{\partial \overline{y}}+\overline{w}\frac{\partial \overline{u}}{\partial \overline{z}}\right]=-{ϵ}^2\frac{\partial \overline{p}}{\partial \overline{x}}+2{ϵ}^2\frac{\partial \overline{\mu }}{\partial \overline{x}}\frac{\partial \overline{u}}{\partial \overline{x}}+\overline{\mu }\left({ϵ}^2\frac{{\partial }^2\overline{u}}{\partial {\overline{x}}^2}+\frac{{\partial }^2\overline{u}}{\partial {\overline{y}}^2}+\frac{{\partial }^2\overline{u}}{\partial {\overline{z}}^2}\right)\end{array}$$

(6a)

$$\begin{array}{c}\hfill \begin{array}{c}\hfill Re{ϵ}^2\left[\frac{\partial \overline{v}}{\partial \overline{t}}+\overline{u}\frac{\partial \overline{v}}{\partial \overline{x}}+\overline{v}\frac{\partial \overline{v}}{\partial \overline{y}}+\overline{w}\frac{\partial \overline{v}}{\partial \overline{z}}\right]=\frac{\partial \overline{\mu }}{\partial \overline{x}}\left(\frac{\partial \overline{u}}{\partial \overline{y}}+{ϵ}^2\frac{\partial \overline{v}}{\partial \overline{x}}\right)-\frac{\partial \overline{p}}{\partial \overline{y}}\\ \hfill +\overline{\mu }\left({ϵ}^2\frac{{\partial }^2\overline{v}}{\partial {\overline{x}}^2}+\frac{{\partial }^2\overline{v}}{\partial {\overline{y}}^2}+\frac{{\partial }^2\overline{v}}{\partial {\overline{z}}^2}\right)\end{array}\end{array}$$

(6b)

$$\begin{array}{c}\hfill \begin{array}{c}\hfill Re{ϵ}^2\left[\frac{\partial \overline{w}}{\partial \overline{t}}+\overline{u}\frac{\partial \overline{w}}{\partial \overline{x}}+\overline{v}\frac{\partial \overline{w}}{\partial \overline{y}}+\overline{w}\frac{\partial \overline{w}}{\partial \overline{z}}\right]=\frac{\partial \overline{\mu }}{\partial \overline{x}}\left(\frac{\partial \overline{u}}{\partial \overline{z}}+{ϵ}^2\frac{\partial \overline{w}}{\partial \overline{x}}\right)-\frac{\partial \overline{p}}{\partial \overline{z}}\\ \hfill +\overline{\mu }\left({ϵ}^2\frac{{\partial }^2\overline{w}}{\partial {\overline{x}}^2}+\frac{{\partial }^2\overline{w}}{\partial {\overline{y}}^2}+\frac{{\partial }^2\overline{w}}{\partial {\overline{z}}^2}\right)\end{array}\end{array}$$

(6c)

while the dynamic and kinematic boundary conditions are transformed into

$$\begin{array}{c}\hfill \begin{array}{c}\hfill {ϵ}^2\left(-\overline{p}+2\overline{\mu }\frac{\partial \overline{u}}{\partial \overline{x}}\right){\overline{n}}_x^{\left(i\right)}+\overline{\mu }\left(\frac{\partial \overline{u}}{\partial \overline{y}}+{ϵ}^2\frac{\partial \overline{v}}{\partial \overline{x}}\right){\overline{n}}_y^{\left(i\right)}+\overline{\mu }\left(\frac{\partial \overline{u}}{\partial \overline{z}}+{ϵ}^2\frac{\partial \overline{w}}{\partial \overline{z}}\right){\overline{n}}_z^{\left(i\right)}=\\ \hfill -{ϵ}^2\left(\overline{\gamma }{\overline{\kappa }}^{\left(i\right)}+{\overline{p}}_H^{\left(i\right)}\right){\overline{n}}_x^{\left(i\right)}\end{array}\end{array}$$

(7a)

$$\left(-\overline{p}+2\overline{\mu }\frac{\partial \overline{v}}{\partial \overline{y}}\right){\overline{n}}_y^{\left(i\right)}+\overline{\mu }\left(\frac{\partial \overline{u}}{\partial \overline{y}}+{ϵ}^2\frac{\partial \overline{v}}{\partial \overline{x}}\right){\overline{n}}_x^{\left(i\right)}+\overline{\mu }\left(\frac{\partial \overline{v}}{\partial \overline{z}}+\frac{\partial \overline{w}}{\partial \overline{y}}\right){\overline{n}}_z^{\left(i\right)}=-\left(\overline{\gamma }{\overline{\kappa }}^{\left(i\right)}+{\overline{p}}_H^{\left(i\right)}\right){\overline{n}}_y^{\left(i\right)}$$

(7b)

$$\left(-\overline{p}+2\overline{\mu }\frac{\partial \overline{w}}{\partial \overline{z}}\right){\overline{n}}_z^{\left(i\right)}+\overline{\mu }\left(\frac{\partial \overline{u}}{\partial \overline{z}}+{ϵ}^2\frac{\partial \overline{w}}{\partial \overline{x}}\right){\overline{n}}_x^{\left(i\right)}+\overline{\mu }\left(\frac{\partial \overline{v}}{\partial \overline{z}}+\frac{\partial \overline{w}}{\partial \overline{y}}\right){\overline{n}}_y^{\left(i\right)}=-\left(\overline{\gamma }{\overline{\kappa }}^{\left(i\right)}+{\overline{p}}_H^{\left(i\right)}\right){\overline{n}}_z^{\left(i\right)}$$

(7c)

and

$$\frac{\partial {\overline{G}}^{\left(i\right)}}{\partial \overline{t}}+\overline{u}\frac{\partial {\overline{G}}^{\left(i\right)}}{\partial \overline{x}}+\overline{v}\frac{\partial {\overline{G}}^{\left(i\right)}}{\partial \overline{y}}+\overline{w}\frac{\partial {\overline{G}}^{\left(i\right)}}{\partial \overline{z}}=0$$

(8)

$ϵ$ and $Re$ are the slenderness parameter and the Reynolds number. They read

$$\begin{array}{c}\hfill ϵ=\frac{{\chi }_0}{L}\end{array}$$

(9a)

$$\begin{array}{c}\hfill Re=\frac{\rho U_{0}L}{{\mu }_s}\end{array}$$

(9b)

2.3. Final Asymptotic Equations

Triggered by the work of previous researchers (Yarin et al. [18], Fitt et al. [19], Luzi et al. [20], Taroni et al. [24], Dewynne et al. [27], Cummings and Howell [30], Stokes et al. [31]), Chen et al. [32] performed a regular expansion of the unknown in terms of the parameter ${ϵ}^2$

$$\begin{array}{cc}\hfill \overline{u}& ={\overline{u}}_0(\overline{t},\overline{x},\overline{y},\overline{z})+{ϵ}^2{\overline{u}}_1(\overline{t},\overline{x},\overline{y},\overline{z})+{ϵ}^4{\overline{u}}_2(\overline{t},\overline{x},\overline{y},\overline{z})\hfill \\ \hfill \overline{v}& ={\overline{v}}_0(\overline{t},\overline{x},\overline{y},\overline{z})+{ϵ}^2{\overline{v}}_1(\overline{t},\overline{x},\overline{y},\overline{z})+{ϵ}^4{\overline{v}}_2(\overline{t},\overline{x},\overline{y},\overline{z})\hfill \\ \hfill \overline{w}& ={\overline{w}}_0(\overline{t},\overline{x},\overline{y},\overline{z})+{ϵ}^2{\overline{w}}_1(\overline{t},\overline{x},\overline{y},\overline{z})+{ϵ}^4{\overline{w}}_2(\overline{t},\overline{x},\overline{y},\overline{z})\hfill \\ \hfill \overline{p}& ={\overline{p}}_0(\overline{t},\overline{x},\overline{y},\overline{z})+{ϵ}^2{\overline{p}}_1(\overline{t},\overline{x},\overline{y},\overline{z})+{ϵ}^4{\overline{p}}_2(\overline{t},\overline{x},\overline{y},\overline{z})\hfill \\ \hfill {\overline{G}}^{\left(i\right)}& ={\overline{G}}_0^{\left(i\right)}(\overline{t},\overline{x},\overline{y},\overline{z})+{ϵ}^2{\overline{G}}_1^{\left(i\right)}(\overline{t},\overline{x},\overline{y},\overline{z})+{ϵ}^4{\overline{G}}_2^{\left(i\right)}(\overline{t},\overline{x},\overline{y},\overline{z})\hfill \end{array}$$

Substituting these expressions into (5), (6a), (6b), (6c), (7a), (7b), (7c), and (8), and assuming that ${ϵ}^{2}Re\ll 1$, Chen et al. [41] obtained at the zero^th order the following mass and momentum equations

$$\frac{\partial {\overline{u}}_0}{\partial \overline{x}}+\frac{\partial {\overline{v}}_0}{\partial \overline{y}}+\frac{\partial {\overline{w}}_0}{\partial \overline{z}}=0$$

(10a)

$$0=\frac{{\partial }^2{\overline{u}}_0}{\partial {\overline{y}}^2}+\frac{{\partial }^2{\overline{u}}_0}{\partial {\overline{z}}^2}$$

(10b)

$$0=\frac{\partial \overline{\mu }}{\partial \overline{x}}\frac{\partial {\overline{u}}_0}{\partial \overline{y}}-\frac{\partial {\overline{p}}_0}{\partial \overline{y}}+\overline{\mu }\left(\frac{{\partial }^2{\overline{v}}_0}{\partial {\overline{y}}^2}+\frac{{\partial }^2{\overline{v}}_0}{\partial {\overline{z}}^2}\right)$$

(10c)

$$0=\frac{\partial \overline{\mu }}{\partial \overline{x}}\frac{\partial {\overline{u}}_0}{\partial \overline{z}}-\frac{\partial {\overline{p}}_0}{\partial \overline{z}}+\overline{\mu }\left(\frac{{\partial }^2{\overline{w}}_0}{\partial {\overline{y}}^2}+\frac{{\partial }^2{\overline{w}}_0}{\partial {\overline{z}}^2}\right)$$

(10d)

From Eq. (10b), the eigensolution ${\overline{u}}_0={\overline{u}}_0\left(\overline{t},\overline{x}\right)$ follows. The kinematic and dynamic boundary conditions at the zero^th order read

$$\frac{\partial {\overline{G}}_0^{\left(i\right)}}{\partial \overline{t}}+{\overline{u}}_0\frac{\partial {\overline{G}}_0^{\left(i\right)}}{\partial \overline{x}}+{\overline{v}}_0\frac{\partial {\overline{G}}_0^{\left(i\right)}}{\partial \overline{y}}+{\overline{w}}_0\frac{\partial {\overline{G}}_0^{\left(i\right)}}{\partial \overline{z}}=0$$

(11a)

$$\frac{\partial {\overline{u}}_0}{\partial \overline{y}}\frac{\partial {\overline{G}}_0^{\left(i\right)}}{\partial \overline{y}}+\frac{\partial {\overline{u}}_0}{\partial \overline{z}}\frac{\partial {\overline{G}}_0^{\left(i\right)}}{\partial \overline{z}}=0$$

(11b)

$$\left(-{\overline{p}}_0+2\overline{\mu }\frac{\partial {\overline{v}}_0}{\partial \overline{y}}\right)\frac{\partial {\overline{G}}_0^{\left(i\right)}}{\partial \overline{y}}+\overline{\mu }\left(\frac{\partial {\overline{v}}_0}{\partial \overline{z}}+\frac{\partial {\overline{w}}_0}{\partial \overline{y}}\right)\frac{\partial {\overline{G}}_0^{\left(i\right)}}{\partial \overline{z}}=-\left(\overline{\gamma }{\overline{\kappa }}_0^{\left(i\right)}+{\overline{p}}_H^{\left(i\right)}\right)\frac{\partial {\overline{G}}_0^{\left(i\right)}}{\partial \overline{y}}$$

(11c)

$$\left(-{\overline{p}}_0+2\overline{\mu }\frac{\partial {\overline{w}}_0}{\partial \overline{z}}\right)\frac{\partial {\overline{G}}_0^{\left(i\right)}}{\partial \overline{z}}+\overline{\mu }\left(\frac{\partial {\overline{v}}_0}{\partial \overline{z}}+\frac{\partial {\overline{w}}_0}{\partial \overline{y}}\right)\frac{\partial {\overline{G}}_0^{\left(i\right)}}{\partial \overline{y}}=-\left(\overline{\gamma }{\overline{\kappa }}_0^{\left(i\right)}+{\overline{p}}_H^{\left(i\right)}\right)\frac{\partial {\overline{G}}_0^{\left(i\right)}}{\partial \overline{z}}$$

(11d)

for $i=0,1,...N$. ${\overline{k}}_0^{\left(i\right)}$ assumes the form

$${\overline{k}}_0^{\left(i\right)}=-\frac{1}{2}\frac{\frac{{\partial }^2{\overline{G}}_0^{\left(i\right)}}{\partial {\overline{y}}^2}+\frac{{\partial }^2{\overline{G}}_0^{\left(i\right)}}{\partial {\overline{z}}^2}}{\sqrt{{\left(\frac{\partial {\overline{G}}_0^{\left(i\right)}}{\partial \overline{y}}\right)}^2+{\left(\frac{\partial {\overline{G}}_0^{\left(i\right)}}{\partial \overline{z}}\right)}^2}}$$

(12)

for more details, see Chen et al. [32]. To obtain the leading-order model for the axial flow, Cummings and Howell [30], Stokes et al. [31], and Chen et al. [32] show that it is necessary to proceed with the asymptotic expansion up to the $\mathcal{O}\left({ϵ}^2\right)$ order

$$Re\left[\frac{\partial {\overline{u}}_0}{\partial \overline{t}}+{\overline{u}}_0\frac{\partial {\overline{u}}_0}{\partial \overline{x}}\right]+\frac{\partial {\overline{p}}_0}{\partial \overline{x}}-2\frac{\partial \overline{\mu }}{\partial \overline{x}}\frac{\partial {\overline{u}}_0}{\partial \overline{x}}-\mu \frac{{\partial }^2{\overline{u}}_0}{\partial {\overline{x}}^2}=\mu \left(\frac{{\partial }^2{\overline{u}}_1}{\partial {\overline{y}}^2}+\frac{{\partial }^2{\overline{u}}_1}{\partial {\overline{z}}^2}\right)$$

(13a)

$$\begin{array}{c}\hfill \overline{\mu }\left(\frac{\partial {\overline{u}}_1}{\partial \overline{y}}\frac{\partial {\overline{G}}_0}{\partial \overline{y}}+\frac{\partial {\overline{u}}_1}{\partial \overline{z}}\frac{\partial {\overline{G}}_0}{\partial \overline{z}}\right)=\left({\overline{p}}_0-2\overline{\mu }\frac{\partial {\overline{u}}_0}{\partial \overline{x}}-\overline{\gamma }{\overline{\kappa }}_0^{\left(i\right)}-{\overline{p}}_H^{\left(i\right)}\right)\frac{\partial {\overline{G}}_0}{\partial \overline{x}}\\ \hfill -\overline{\mu }\left(\frac{\partial {\overline{v}}_0}{\partial \overline{x}}\frac{\partial {\overline{G}}_0}{\partial \overline{y}}+\frac{\partial {\overline{w}}_0}{\partial \overline{x}}\frac{\partial {\overline{G}}_0}{\partial \overline{z}}\right)\end{array}$$

(13b)

Integrating now Eqs. (13a) and (13b) in the cross-sectional plane and applying the divergence theorem, one obtains after lengthy algebraic manipulations [31,32]

$$\frac{\partial }{\partial \overline{x}}\left(3\overline{\mu }\frac{\partial u_0}{\partial x}{\overline{S}}_0\right)=Re{\overline{S}}_0\left[\frac{\partial {\overline{u}}_0}{\partial \overline{t}}+{\overline{u}}_0\frac{\partial {\overline{u}}_0}{\partial \overline{x}}\right]-\frac{1}{2}\overline{\gamma }\frac{\partial {\overline{\mathrm{\Gamma }}}_0}{\partial x}+\frac{\partial {\overline{p}}_H^{\left(0\right)}}{\partial \overline{x}}{\overline{S}}_0^0-\sum _{k=1}^N\frac{\partial {\overline{p}}_H^{\left(i\right)}}{\partial \overline{x}}{\overline{S}}_0^{\left(k\right)}$$

(14)

where

$${\overline{S}}_0={\overline{S}}_0^0-\sum _{i=1}^N{\overline{S}}_0^{\left(i\right)}$$

is the difference between the area encompassed by the external boundary and the sum of the areas encompassed by the internal boundaries. Since $Re\ll 1$ the inertia term can be neglected. In addition, we consider the case where the internal channel pressurization is constant. Therefore, Eq. (14) reduces to

$$\frac{\partial }{\partial \overline{x}}\left(3\overline{\mu }\frac{\partial u_0}{\partial x}{\overline{S}}_0\right)=-\frac{1}{2}\overline{\gamma }\frac{\partial {\overline{\mathrm{\Gamma }}}_0}{\partial x}$$

(15)

Eq. (15) can be written in a Lagrangian reference frame where the problem is unsteady as

$$3\overline{\mu }\frac{D{\overline{S}}_0}{D\overline{t}}-\frac{1}{2}\overline{\gamma }{\overline{\mathrm{\Gamma }}}_0=-6\overline{\sigma }$$

(16)

where $D{\overline{S}}_0/D\overline{t}=\partial {\overline{S}}_0/\partial \overline{t}+{\overline{u}}_0\partial {\overline{S}}_0/\partial \overline{x}$ is the material derivative and $6\overline{\sigma }$ is the scaled pulling tension of the fibre $6\overline{\sigma }=\sigma L/\left({\mu }_0{U}_0{S}_0\right)$. To this end, Stokes et al. [31] integrated the kinematic condition Eq. (11a) in the cross-sectional plane and utilized the zero^th order continuity equation,i.e., Eq. (10a).

2.3.1. Leading-order model for the transverse flow

Stokes et al. [31], Stokes et al. [33] and Chen et al. [32] proposed the following coordinate system

$$\begin{array}{cccccccccccc}\hfill \overline{t}& =\tilde{t},\hfill & \hfill \overline{x}& ={\int }_0^{\tilde{t}}{\overline{u}}_0\left(\tilde{x},T\right)dT+\tilde{x},\hfill & \hfill \overline{y}& =\overline{\chi }\tilde{y},\hfill & \hfill \overline{z}& =\overline{\chi }\tilde{z},\hfill & \hfill \overline{\mathrm{\Gamma }}& =\overline{\chi }\tilde{\mathrm{\Gamma }},\hfill & {\tilde{\kappa }}^{\left(i\right)}& =\overline{\chi }{\overline{\kappa }}^{\left(i\right)}\end{array}$$

(17)

for the transverse flow in the cross-sectional plane. By doing so, neither rigid-body translation nor rotation of each cross-section is assumed. The pressure and velocity components in the cross-sectional plane in this new coordinate system transform as follows

$$\begin{array}{cccccc}\hfill {\overline{p}}_0& =-\overline{\mu }\frac{\partial {\overline{u}}_0}{\partial \overline{x}}+\frac{\overline{\gamma }}{\overline{\chi }}\tilde{p},\hfill & \hfill {\overline{v}}_0& =-\frac{1}{2}\frac{\partial {\overline{u}}_0}{\partial \overline{x}}\overline{y}+\frac{\overline{\gamma }}{\overline{\mu }}\tilde{v},\hfill & \hfill {\overline{w}}_0& =-\frac{1}{2}\frac{\partial {\overline{u}}_0}{\partial \overline{x}}\overline{z}+\frac{\overline{\gamma }}{\overline{\mu }}\tilde{w}\hfill \end{array}$$

(18)

In turn, the leading-order equations for the mass, momentum, kinematic, and dynamic boundary conditions read

$$\frac{\partial \tilde{v}}{\partial \tilde{y}}+\frac{\partial \tilde{w}}{\partial \tilde{z}}=0$$

(19a)

$$\frac{{\partial }^2\tilde{v}}{\partial {\tilde{y}}^2}+\frac{{\partial }^2\tilde{v}}{\partial {\tilde{z}}^2}=\frac{\partial \tilde{p}}{\partial \tilde{y}}$$

(19b)

$$\frac{{\partial }^2\tilde{w}}{\partial {\tilde{y}}^2}+\frac{{\partial }^2\tilde{w}}{\partial {\tilde{z}}^2}=\frac{\partial \tilde{p}}{\partial \tilde{z}}$$

(19c)

$$\frac{\partial {\overline{G}}^{\left(i\right)}}{\partial \overline{\tau }}+\tilde{v}\frac{\partial {\overline{G}}^{\left(i\right)}}{\partial \tilde{y}}+\tilde{w}\frac{\partial {\overline{G}}^{\left(i\right)}}{\partial \tilde{z}}=0$$

(19d)

$$\left(-\tilde{p}+2\frac{\partial \tilde{v}}{\partial \tilde{y}}\right)\frac{\partial {\tilde{G}}^{\left(i\right)}}{\partial \tilde{y}}+\left(\frac{\partial \tilde{v}}{\partial \tilde{z}}+\frac{\partial \tilde{w}}{\partial \tilde{y}}\right)\frac{\partial {\tilde{G}}^{\left(i\right)}}{\partial \tilde{z}}=-\left({\tilde{\kappa }}^{\left(i\right)}+\frac{{\overline{p}}_H^{\left(i\right)}}{\overline{\gamma }}\overline{\chi }\right)\frac{\partial {\tilde{G}}^{\left(i\right)}}{\partial \tilde{y}}$$

(19e)

$$\left(-\tilde{p}+2\frac{\partial \tilde{w}}{\partial \tilde{z}}\right)\frac{\partial {\tilde{G}}^{\left(i\right)}}{\partial \tilde{z}}+\left(\frac{\partial \tilde{v}}{\partial \tilde{z}}+\frac{\partial \tilde{w}}{\partial \tilde{y}}\right)\frac{\partial {\tilde{G}}^{\left(i\right)}}{\partial \tilde{y}}=-\left({\tilde{\kappa }}^{\left(i\right)}+\frac{{\overline{p}}_H^{\left(i\right)}}{\overline{\gamma }}\overline{\chi }\right)\frac{\partial {\tilde{G}}^{\left(i\right)}}{\partial \tilde{z}}$$

(19f)

where Eqs. (19d), (19e), and (19f) are defined on ${\tilde{G}}^{\left(i\right)}(\overline{\tau },\tilde{y},\tilde{z})=0$. (19d) is written in terms of a “reduced time” $\overline{\tau }$ defined as

$$\overline{\tau }=\overline{\gamma }{\int }_0^{\tilde{t}}\frac{d\overline{T}}{\overline{\mu }\overline{\chi }}$$

(20)

which is referred to as the total deformation time of the draw [30,31,32]. Eqs. (19d), (19e), and (19f) constitutes a Stokes-flow free-boundary problem defined in a region with unit viscosity, surface tension, and cross-sectional area. Using Eq. (20), Eq. (16) becomes

$$\frac{\partial \overline{\chi }}{\partial \overline{\tau }}-\frac{\overline{\chi }}{12}\tilde{\mathrm{\Gamma }}\left(\overline{\tau }\right)=-\frac{\overline{\sigma }}{\overline{\gamma }}$$

(21)

for more details see [30,31,32,33].

2.3.2. The case of a circular tube

The asymptotic mass and momentum equations are first employed to model the drawing of capillaries. Chen et al. [32] naturally chose polar coordinates $(r,\theta )$ to recast Eqs. (19a), (19b), (19c), (19d), (19e), and (19f) for the cross-sectional plane flow. The outer and inner radii are indicated by $\tilde{r}=\tilde{R}\left(\overline{\tau }\right)$ and $\tilde{r}=\tilde{\rho }\left(\overline{\tau }\right)\tilde{R}\left(\overline{\tau }\right)$, respectively. $\tilde{\rho }\left(\overline{\tau }\right)$ is the ratio between the inner and the outer radii such that

$$\pi {\tilde{R}}^2\left(1-{\tilde{\rho }}^2\right)=1$$

The axial momentum equation, Eq. (21), transforms as follows

$$\frac{\partial \overline{\chi }}{\partial \overline{\tau }}=\frac{\sqrt{\pi }}{6}\overline{\chi }\sqrt{\frac{1+\tilde{\rho }}{1-\tilde{\rho }}}-\frac{\overline{\sigma }}{\overline{\gamma }}$$

(22)

and the cross-plane equations read

$$\begin{array}{ccccc}\hfill \frac{1}{\tilde{r}}\frac{\partial }{\partial \tilde{r}}\left(\tilde{r}\frac{\partial {\tilde{v}}_r}{\partial \tilde{r}}\right)& =0,\hfill & \hfill \mathrm{in}& \hfill & \hfill \overline{S},\end{array}$$

(23a)

$$\begin{array}{cccccc}\hfill -\frac{\partial \tilde{p}}{\partial \tilde{r}}+\frac{1}{\tilde{r}}\frac{\partial }{\partial \tilde{r}}\left(\tilde{r}\frac{\partial {\tilde{v}}_r}{\partial \tilde{r}}\right)-\frac{{\tilde{v}}_r}{{\tilde{r}}^2}& =0,\hfill & \hfill & \mathrm{in}\hfill & \hfill & \overline{S},\hfill \end{array}$$

(23b)

$$\begin{array}{cccccc}\hfill 2\frac{\partial {\tilde{v}}_r}{\partial \tilde{r}}-\tilde{p}& =-\left(\tilde{k}+\frac{\overline{\chi }}{\overline{\gamma }}{\overline{p}}_H^{\left(1\right)}\right),\hfill & \hfill & \mathrm{on}\hfill & \hfill & {\tilde{G}}^{\left(1\right)}=0,\hfill \end{array}$$

(23c)

$$\begin{array}{cccccc}\hfill 2\frac{\partial {\tilde{v}}_r}{\partial \tilde{r}}-\tilde{p}& =-\tilde{k},\hfill & \hfill & \mathrm{on}\hfill & \hfill & {\tilde{G}}^{\left(0\right)}=0,\hfill \end{array}$$

(23d)

$$\begin{array}{cccccccccc}\hfill \frac{\partial {\tilde{G}}^{\left(i\right)}}{\partial \overline{\tau }}+{\tilde{v}}_r& =0,\hfill & \hfill & \mathrm{on}\hfill & \hfill & {\tilde{G}}^{\left(i\right)}=0,\hfill & \hfill & \mathrm{for}\hfill & \hfill & i=0,1,\hfill \end{array}$$

(23e)

where $\tilde{G}={\tilde{G}}_O+{\tilde{G}}_I$, for more details see [31,32]. This system of equations can be solved for the radial component of the velocity and the pressure, and an evolution equation for the ratio between the inner and the outer diameter can be found

$${\tilde{v}}_r=\frac{\tilde{\rho }\tilde{R}}{2\tilde{r}\left(\tilde{\rho }-1\right)}\left(1-\frac{{\overline{p}}_H^{\left(1\right)}}{\overline{\gamma }}\frac{\tilde{\rho }\tilde{R}\overline{\chi }}{\tilde{\rho }+1}\right)$$

(24a)

$$\tilde{p}=\frac{1}{\left(1-\tilde{\rho }\right)\tilde{R}}-\frac{\overline{\chi }}{\overline{\gamma }}\frac{{\tilde{\rho }}^2}{\left(1-{\tilde{\rho }}^2\right)}{\overline{p}}_H^{\left(i\right)}$$

(24b)

$$\frac{d\tilde{\rho }}{d\tilde{\tau }}=-\frac{\sqrt{\pi }}{2}{\left(1+\tilde{\rho }\right)}^{3/2}{\left(1-\tilde{\rho }\right)}^{1/2}+\frac{1}{2}\frac{{\overline{p}}_H^{\left(1\right)}}{\overline{\gamma }}\tilde{\rho }\overline{\chi }$$

(24c)

for more details, see Chen et al. [32].

2.3.3. The cross-plane problem and the complex variable formulation

The Stokes equations govern the steady free-boundary cross-plane problem

$$-\tilde{\nabla }\tilde{p}+\overline{\mu }{\tilde{\nabla }}^2\tilde{\mathit{v}}=0$$

(25a)

$$\tilde{\nabla }·\tilde{\mathit{v}}=0$$

(25b)

where the kinematic boundary conditions read

$$\tilde{\mathit{v}}·{\tilde{\mathit{n}}}^{\left(i\right)}=V_n^{\left(i\right)}$$

(25c)

and the dynamic boundary conditions assume the form

$$-\tilde{p}{\tilde{\mathit{n}}}^{\left(i\right)}+2\overline{\mu }\tilde{\mathit{e}}·{\tilde{\mathit{n}}}^{\left(i\right)}=-\overline{\gamma }\tilde{k}{\tilde{\mathit{n}}}^{\left(i\right)}-{\overline{p}}_H^{\left(i\right)}{\tilde{\mathit{n}}}^{\left(i\right)}$$

(25d)

where $\tilde{\mathit{v}}$, $\tilde{\nabla }$, and ${\tilde{\nabla }}^2$, are the velocity vector, gradient, and Laplacian in the two-dimensional cross-sectional plane. $V_n^{\left(i\right)}$ is the normal velocity of the interface and $\tilde{\mathit{e}}$ the cross-plane stress tensor [43]. It can be shown that the two-dimensional problem described by Eqs. (25a) and (25b) can be formulated in terms of a stream function $\psi (y,z)$ that satisfies

$${\tilde{\nabla }}^4\psi =0$$

whose solution is

$$\psi =Im\left[\overline{s}f(s,\overline{\tau })+g(s,\overline{\tau })\right]$$

being $f(s,\overline{\tau })$ and $g(s,\overline{\tau })$ two Goursat functions that are analytic in the fluid region, and $s=\tilde{y}+i\tilde{z}$ and $\overline{s}=\tilde{y}-i\tilde{z}$ are the standard complex variables. The fluid pressure $\tilde{p}$, the vorticity $\tilde{\omega }$, the velocity $\tilde{\mathit{v}}$, and the rate-of-strain tensor $\tilde{\mathit{e}}$ may be written as

$$\frac{\tilde{p}}{\overline{\mu }}-i\tilde{\omega }=4f^{\prime }\left(s,\overline{\tau }\right)$$

(26a)

$$\tilde{v}+i\tilde{w}=-f\left(s,\overline{\tau }\right)+s\overline{f^{\prime }(s,\overline{\tau })}+\overline{g^{\prime }(s,\overline{\tau })}$$

(26b)

$$e_{11}+e_{12}=s\overline{f^{\prime \prime }(s,\overline{\tau })}+\overline{g^{\prime \prime }(s,\overline{\tau })}$$

(26c)

The dynamic boundary condition Eq. (25d) assume the form

$$f\left(s,\overline{\tau }\right)+s\overline{f^{\prime }(s,\overline{\tau })}+\overline{g^{\prime }(s,\overline{\tau })}=-\frac{1}{2}i\frac{\overline{\gamma }}{\overline{\mu }}\frac{ds}{d\theta }+\frac{1}{2}s\frac{{\overline{p}}_H^{\left(k\right)}}{\overline{\mu }}$$

(26d)

for $k=1,...,N$, and the kinematic boundary condition reads

$$Re\left\{\frac{ds}{d\overline{\tau }}i\frac{d\overline{s}}{d\theta }\right\}=Re\left\{-f\left(s,\overline{\tau }\right)+s\overline{f^{\prime }(s,\overline{\tau })}+\overline{g^{\prime }(s,\overline{\tau })}i\frac{d\overline{s}}{d\theta }\right\}$$

(26e)

where $\theta$ is the arclength of the boundary $s\left(\theta \right)$.

2.3.4. The Generalised Elliptical Pore Model (GEPM)

Buchak et al. [34] proposed the Generalized Elliptical Pore Model (GEPM) to simulate the drawing process of multi-hole fibres. It assumes that the n channels of a MOF are always elliptical during the evolution and exert a pressure $p_n$, vorticity ${\omega }_n$, and strain rate $k_n$ on the flow field. The Goursat functions $f(s,\overline{\tau })$ and $g(s,\overline{\tau })$ are defined as

$$f\left(s,\overline{\tau }\right)=\sum _{n=1}^M\frac{{\lambda }_n\left(\overline{\tau }\right)}{2\pi (s-L_n\left(\overline{\tau }\right))}+\frac{1}{4}\tilde{P}\left(\overline{\tau }\right)s$$

(27a)

and

$$g^{\prime }\left(s,\overline{\tau }\right)=\sum _{n=1}^M\frac{m_n\left(\overline{\tau }\right)}{2\pi (s-L_n\left(\overline{\tau }\right))}+\sum _{n=1}^M\frac{{\lambda }_n\left(\overline{\tau }\right)\overline{L_n\left(\overline{\tau }\right)}}{2\pi {(s-L_n\left(\overline{\tau }\right))}^2}$$

(27b)

so that the whole fluid flow region is described by a superposition of sources of strength $m_n\left(\overline{\tau }\right)$ and stresslets of strength ${\lambda }_n\left(\overline{\tau }\right)$ that model the presence of the channels. In turn, the pressure term of Eq. (27a) describes the influence of the outer free surface on the inner channels. In Eqs. (27a) and (27b) $L_n$ are the centroids of the nth channels. Buchak et al. [34] derived evolution equations for the channel centroids

$$\overline{\frac{d{L}_n}{d\overline{\tau }}}=\overline{u\left(L_n\right)}=\sum _{j\ne n}^M\frac{m_j}{2\pi (L_n-L_j)}-\sum _{j\ne n}^M\frac{{\lambda }_j\overline{(L_n-L_j)}}{2\pi {(L_n-L_j)}^2}-\overline{{\sum }_{j\ne n}^M\frac{{\lambda }_j}{2\pi (L_n-L_j)}}$$

(28)

and showed that the local pressure, vorticity, and strain rate exerted by the nth channel obey the following relationships

$$\frac{1}{4}\left({\tilde{p}}_n-i{\tilde{\omega }}_n\right)=\frac{1}{4}\tilde{P}-\sum _{j\ne n}^M\frac{{\lambda }_j}{2\pi {(L_n-L_j)}^2}$$

(29a)

$$k_n=-\sum _{j\ne n}^M\frac{m_j}{2\pi {(L_n-L_j)}^2}-\sum _{j\ne n}^M\frac{{\lambda }_j\overline{(L_j-L_n)}}{\pi {(L_n-L_j)}^3}$$

(29b)

Assuming an elliptical shape of the channels in terms of complex variables

$$s_n(\zeta ,\overline{\tau })=L_n\left(\overline{\tau }\right)+\frac{{\alpha }_n\left(\overline{\tau }\right)}{\zeta }+{\beta }_n\left(\overline{\tau }\right)\zeta$$

(30)

Buchak et al. [34] obtained evolution equations for the parameters ${\alpha }_n\left(\overline{\tau }\right)$ and ${\beta }_n\left(\overline{\tau }\right)$ utilizing the method of inner-outer matched asymptotic expansion

$$\frac{\partial {\alpha }_n}{\partial \overline{\tau }}=-{\alpha }_n{I}_n\left(0\right)-\frac{1}{2}{\alpha }_n\left({\tilde{p}}_n-{\overline{p}}_H^{\left(n\right)}\frac{\overline{\chi }}{\overline{\gamma }}\right)$$

(31a)

$$\frac{\partial {\beta }_n}{\partial \overline{\tau }}=-{\beta }_n{I}_n\left(0\right)+\frac{1}{2}{\beta }_n\left({\tilde{p}}_n-{\overline{p}}_H^{\left(n\right)}\frac{\overline{\chi }}{\overline{\gamma }}\right)+2{\overline{k}}_n{\alpha }_n+i{\tilde{\omega }}_n{\beta }_n$$

(31b)

where

$$\frac{m_n}{\pi }=\frac{d}{d\overline{\tau }}\left({\alpha }_n^2-\left|{\beta }_n^2\right|\right)$$

(31c)

$$\frac{{\lambda }_n}{\pi }={\beta }_n\frac{\partial {\alpha }_n}{\partial \overline{\tau }}-{\alpha }_n\frac{\partial {\beta }_n}{\partial \overline{\tau }}+i{\alpha }_n{\beta }_n\widetilde{{\omega }_n}+2{\alpha }_n{\beta }_n{I}_n\left(0\right)-\frac{1}{2}{\alpha }_n^2{I}_n^{{ }^{\prime \prime }}\left(0\right)$$

(31d)

where $m_n\in \mathbb{R}$ and ${\lambda }_n\in \mathbb{C}$. ${\alpha }_n\left(\overline{\tau }\right)\in \mathbb{R}$ and ${\beta }_n\left(\overline{\tau }\right)\in \mathbb{C}$ describe the orientation, eccentricity, and area of each channel during the evolution. The outer boundary is modelled as $\left|s\right|=\tilde{R}\left(\overline{\tau }\right)$ and stays circular and centred at the origin during the evolution. Its evolution equation is

$$\frac{d\tilde{R}}{d\overline{\tau }}=\frac{M}{2\pi \tilde{R}}$$

(31e)

where $M={\sum }_{j=1}^n{m}_n$, for more details see Buchak et al. [34].

2.3.5. Fibre temperature profile and glass viscosity

We select the following functional form of the fibre temperature distribution

$$\overline{T}\left(\overline{x}\right)={\overline{T}}_M\left(a+b*exp\left(c*{\left(\overline{x}-0.5\right)}^2\right)\right)$$

(32)

where ${\overline{T}}_M=T_M/T_S$. $T_M$ is the furnace peak temperature and $T_S=1900K$ represents the softening temperature of the glass. Besides, $a=0.5$, $b=0.4454$, $c=-3.032$, and $d=0.497$. By doing so, the fibre temperature distribution is centred in the middle of the “hot zone”, that is a small region of the oven where the temperature is high enough so that the glass melts and the drawing can occur. Moreover, we retain the length of the “hot zone” to be 12 $\left[mm\right]$, as we did in our previous work. An example of the fibre temperature profile is presented in Figure 2, where the oven peak temperature is $T_M=$ 2050 ${[}^{\circ }C]$.

To set the dependence of the viscosity on the glass temperature, we select

$$\mu =5.8·{10}^{-8}exp\left(\frac{515400}{8.3145T}\right)$$

(33a)

if the temperature range is ${1400}^{\circ }C\le T\le {2500}^{\circ }C$. Instead, we utilize

$$\mu =3.8·{10}^{-14}exp\left(\frac{712000}{8.3145T}\right)$$

(33b)

if the temperature range is ${1000}^{\circ }C\le T\le {1400}^{\circ }C$. In Eqs. (33a) and (33b), the temperature T is in K and the viscosity $\mu$ is in $Pa·s$ as explained by Voyce et al. [23].

3. Results

3.1. Solution methodology

We employ the robust fourth-order Runge-Kutta-Merson method to integrate Eq. (24c) for the case of a single capillary and Eqs. (28), (31a), and (31b) for the case of multi-hole fibres. The parameters ${\tilde{p}}_n$, ${\tilde{\omega }}_n$, $k_n$, $m_n$, and ${\lambda }_n$ are found by solving the algebraic equations Eqs. (29a), (29b), (31c), and (31d). The ratios DR 36-1, DR 54-15, and DR 72-2 are termed slow drawing ratios, and the ratios DR 1-102, DR 2-204, DR 3-306, DR 20-187, and DR 10-42 are termed high drawing ratios. The main initial values and the parameters used for the numerical computations are listed in Table 1.

The values of the glass parameters are in principle all dependent on the temperature. However, we only consider the dependency of the glass viscosity on the temperature since it is the only one that spans several orders of magnitude. In this work, all the governing differential equations are integrated with respect to $\overline{x}$ and not with respect to $\overline{\tau }$, using a step size $\Delta \overline{x}=2.56e^{-4}$. By doing so, the simulations are very stable and it is not necessary to integrate in a loop until the final stage of the drawing is achieved [33]. Instead, we integrate the system of equations iteratively until the value of the tension satisfies the condition

$$\frac{{\overline{u}}_0(\overline{x}=1)}{{\overline{u}}_F}-1<{10}^{-3}$$

3.2. Annular Capillaries: Slow Drawing Ratios (SDRs)

We first consider the case DR 36-1 without internal pressurization to show the evolution of the axial velocity ${\overline{u}}_0$ and the inner and outer fibre radii ${\overline{h}}_1$ and ${\overline{h}}_2$ in the drawing direction $\overline{x}$. Figure 3 a) shows that ${\overline{u}}_0$ varies by several orders of magnitude in the “hot-zone”. At the beginning of the“hot-zone”, the values of the axial velocity ${\overline{u}}_0$ are very low and close to the ones of the feeding speed $U_F$ until approximately $\overline{x}$≈ 0.2. Subsequently, ${\overline{u}}_0$ increases until it achieves the final value at approximately $\overline{x}$≈ 0.9. The main variations of the values of the axial velocity occur in the part of the “hot-zone” where the viscosity is low enough to allow the molten glass to flow. To obey the mass conservation equation, the fibre-cross section reduces its size, from values very close to $\overline{\chi }$≈ 1 at the beginning of the “hot-zone” down to $\overline{\chi }$≈ 0.06 at the end of it. In Figure 4 we compare the external diameters of an annular capillary obtained numerically and experimentally for the three drawing ratios SDR 36-1, SDR 54-15, and SDR 72-2 without internal pressurization. We term these drawing ratios “slow drawing ratios” (SDRs), since the drawing speeds are significantly lower compared to the ones needed to draw fibres of commercially available size. We progressively increment the oven peak temperature from 1850 ${[}^{\circ }C]$ up to 2050 ${[}^{\circ }C]$ in steps of 25 ${[}^{\circ }C]$. The differences between numerical and experimental results are very small since the discrepancies vary from a minimum of approximately 0.06 % up to a maximum of approximately 3.5 %, see Figure 4 a), b), and c). In Figure 5 we contrast the air-filling fractions obtained numerically and experimentally for the same SDRs considered in 4. In this case, the numerical results only agree reasonably well with experimental data in the peak temperature range between 1850 ${[}^{\circ }C]$ and 2000 ${[}^{\circ }C]$, since the differences between them vary approximately from 0.14 % up to a maximum of approximately 17.5 %. In the case of higher oven peak temperatures, the discrepancies between numerical and experimental data can reach very high values. For instance, in the case of DR 36-1 and $T_{peak}=$ 2050 ${[}^{\circ }C]$ the difference between the air-filling fraction obtained numerically and experimentally is approximately 61.2 %. Since the fibre temperature correlation we propose in this work must encompass a wide range of peak temperatures, the fibre temperature is probably underestimated in the late stage of the drawing, especially in the case of high furnace peak temperatures. This leads to an overprediction of the viscosity of the glass and consequently, of the fibre size.

Afterward, in Figure 6 a) and Figure 6 b), we analyze a case with internal pressurization. We fix the furnace peak temperature $T_{peak}=$ 1950 ${[}^{\circ }C]$ and vary the internal pressure $p_H$ from 5 $\left[mbar\right]$ to 25 $\left[mbar\right]$ in steps of 2.5 $\left[mbar\right]$. In Figure 6a) we contrast the final external diameter and in Figure 6b) the air-filling fractions of an annular capillary obtained both experimentally and numerically. As far as the external diameter of the capillary is concerned, the numerical data agree well with the experimental ones for the whole range of pressure values considered, since the maximum difference between experiments and computations is approximately 12 %. With regard to the air-filling fractions, the concordance between experiments and simulations is very good in the pressure range 5 $\left[mbar\right]$$<p_H<$ 17.5 $\left[mbar\right]$, since the discrepancies between numerical and experimental values achieve a maximum of approximately 18 %. In the pressure range 20 $\left[mbar\right]$$<p_H<$ 25 $\left[mbar\right]$ the agreement between experiments and simulations declines but it is still fair, reaching a maximum value of approximately 24 %. The major deviations between experimental and numerical results occur as the value of the inner pressurization increases. This may be caused by the fibre temperature functional form used in this work, which may lead to an overestimation of the viscosity of the glass for the considered oven peak temperature. Consequently, the inner surface of the capillary computed numerically results less inflated than the real one under the effects of internal pressurization, and the final fibre size results under-predicted.

3.3. Annular Capillaries: High Drawing Ratios (HDRs)

In this section, we analyze the cases of three high drawing ratios (HDS), that is, DR 1-102, DR 2-204, and DR 3-306, both in the presence and in the absence of internal pressurization. In these cases, the final size of an annular capillary reaches dimensions similar to those of commercially available standard optical fibres. The evolution of the axial component of the velocity, of the fibre cross-section, and of the inner and outer radii are similar to those presented in Figure 3. Nevertheless, the variation between the initial and the final values is much steeper and occurs in a shorter distance [25]. As we did previously, we begin the comparison between experiments and simulations considering two cases without internal pressurization, that is, the drawing ratios DR 1-102 and DR 3-306. To this end, we progressively change the furnace peak temperature from $T_{peak}=$ 1900 ${[}^{\circ }C]$ to $T_{peak}=$ 2050 ${[}^{\circ }C]$. The agreement between the final outer diameter of the capillary obtained numerically and experimentally is very good, and the maximum difference between the experimental and numerical results is approximately 12 %. This happens for the case DR 1-102 and $T_{peak}=$ 2050 ${[}^{\circ }C]$, where the capillary may have undergone experimentally a partial collapse, see Figure 7 a). With respect to the air-filling fraction, the agreement between experiments and numerical simulations is also very good for all the cases of the oven peak temperature considered, except for the case DR 3-306 and $T_{peak}=$ 2050 ${[}^{\circ }C]$, where the deviation is approximately 24.5 %, see Figure 8 a). Analogously to the SDRs cases without internal pressurization, the deviations between experiments and simulations increase as the furnace peak temperature increases, and the functional form of the fibre temperature profile we selected may under-predict the real one. Therefore, the predicted final fibre size is smaller than the real one due to the excessive values the viscosity assumes.

In Figure 7 c) and Figure 8 c), we consider a case with internal pressurization. Specifically, we define a case with DR 2-204, $T_{peak}=$ 1950 ${[}^{\circ }C]$, and inner pressurization ranging from $p_0=5\left[mbar\right]$ up to $p_0=40\left[mbar\right]$. The agreement between experiments and numerical results are very good for the whole pressure range, and the maximum discrepancy between the numerical and experimental outer diameter is approximately 13.8 %, see Figure 7 c). The agreement between the numerical and experimental data of the air-filling fractions is also very good and the maximum difference is approximately 15.5 %, see Figure 8 c).

3.4. Holey Fibres (HFs)

In this section, we first test our fibre temperature profile for the case of a six-hole structure in the absence of internal pressurization. We previously simulated such a geometrical configuration utilizing Ansys PolyFlow^®, a FEM-based commercial software. The initial hole sizes and configuration can be found in Fig. 1 of the manuscript of Luzi et al. [40]. The drawing conditions are: $T_{peak}=$ 1890 ${[}^{\circ }C]$, DR 20-187, and $p_H=$ 0 $\left[mbar\right]$. In Figure 9a), we overlap the experimental results, the previous FEM simulations, and the outcome of the present computations. Figure 9a) exhibits an excellent agreement between the current computations and the experimental results. The shape, size, and positions of the holes at the end of the drawing are well predicted, and the maximum discrepancy between simulations and experiments is only 0.17 %. Afterward, we consider a thirty-six-hole cross-section with internal pressurization, whose geometrical details of the shape are given in the manuscript of Frosz et al. [44]. Besides, $T_{peak}=$ 1973 ${[}^{\circ }C]$, $p_H=$ 173 $\left[mbar\right]$, and DR 10-42 as indicated in Table 1. Figure 9b) displays an excellent agreement between the final fibre cross-section obtained experimentally and the one computed numerically. In particular, the hole sizes and locations in the fibre cross-section are correctly predicted by the numerical simulations.

4. Discussion

We proposed a novel correlation to compute the fibre temperature profile and integrated it into the three-dimensional asymptotic extensional-flow models of Stokes et al. [31], Chen et al. [32] and Buchak et al. [34] assuming thermal equilibrium between the furnace walls and the fibre. Those models describe the drawing process of both annular capillaries and MOFs that can contain many holes. In our previous work [25], we showed that the temperature profile of a capillary of certain size strongly varies in the drawing direction and imperceptible changes can be noted in the radial one. Besides, in the first half of the “hot zone" the fibre and furnace temperature profile almost coincided both for the SDRs cases and for the HDRs cases, indicating thermal equilibrium between the oven and the capillary due to the low values of the feeding velocities. In the second half of the “hot zone" the fibre temperature profile significantly departed from that of the oven, assuming always lower values, thereby revealing a strong convective cooling phenomenon due to the high values of the drawing speeds. This drove us toward the selection of a fibre temperature profile with a peak temperature lower compared to the one of the furnace but with a slightly wider width to avoid excessive high values of the viscosity that could prejudice the drawing leading to the breakage of the fibre. After a trial and error procedure, we fixed the values of the parameters that determine the fibre temperature profile. We chose to integrate the system of differential equations with respect to $\overline{x}$ and not with respect to $\overline{\tau }$. By doing so, all the numerical simulations are very stable and there is no need to implement an upper bound for the viscosity in the computer codes.

We first compared the results of our numerical simulations for the case of annular capillaries and afterward we considered holey fibres. To this end, we contrasted our numerical outcomes with the experiments of Luzi et al. [20], Luzi et al. [40] and Frosz et al. [44]. Overall, we obtained good agreements between numerical and experimental results for a wide range of drawing conditions and fibre cross-sectional shapes. Besides, the computational time needed to solve the governing system of differential and algebraic equations is substantially shorter than the time needed to solve the full set of three-dimensional Navier-Stokes equations. On the one hand, it takes some minutes to solve numerically the equations governing the drawing of annular capillaries and to perform several iterations for the case with thirty-six capillaries. On the other hand, it takes hours on a High Performance Computer (HPC) center to carry out a simulation of a six-hole fibre structure using a FEM-based sofwtare [40]. The fibre temperature profile utilized in this work leads to very good agreement between numerical simulations and experimental data for the cases of HDRs and furnace peak temperatures $T_{peak}\le$ 2000 ${[}^{\circ }C]$. Plausible reasons for this are the lower peak temperature of the fibre compared to the one of the oven and the effects of the strong convective cooling. Both may lead to values of the viscosity that are so high that the drawing conditions only influence to a certain extent the final fibre cross-section. Differently, for the cases of SDRs and furnace peak temperatures $T_{peak}\ge$ 2000 ${[}^{\circ }C]$ the agreement between numerical computations and experiments depreciate significantly. The values of the viscosity may be low enough to allow significant changes of the fibre geometry during the drawing process. Besides, the convective cooling effects may be insufficient due to the low drawing speeds. In addition, a partial collapse of the fibre structure begin to occur in the experiments due to longer residence time of the fibre in the oven and a longer exposition of the fibre to the furnace heating, modifying significantly the fibre geometry. Under those circumstances, the proposed fibre temperature profile does not lead to the correct prediction of the glass viscosity which may be low enough to enable significant changes of the fibre structure during the drawing. Consequently, the final fibre size result overestimated. The specified fibre temperature profile also lead to values of the final fibre structure that slightly change with the oven peak temperature. We hypothesize that the values of the viscosity computed using the proposed correlation for the temperature of the fibre are too high to allow the fibre to change significantly its size with the oven peak temperature.

To obtain more accurate results for a wide range of furnace peak temperatures and drawing conditions the heat exchange between oven and furnace must be modelled [24,25,33]. The fibre temperature profile strongly depends on the heat exchange between the oven and the furnace and on the convective cooling, in particular in the late stage of the drawing where the fibre achieves high speeds. Besides, the effects of viscous dissipation may also become important in regions where the velocity gradients are significant, especially in the axial direction, and capillaries may begin to collapse at high peak temperatures if the feeding and drawing speeds are low enough. In such circumstances, first the inner surface of a capillary begins to collapse and afterward the whole fibre structure. All these issues make the final fibre structure difficult to accurately predict with only one correlation that solely depends on the peak temperature of the oven. Nevertheless, we obtain a good concordance between experiments and simulations also for the case of holey fibres, both when the fibre is not pressurized and when internal pressurization is applied. Since in this case the furnace the peak temperature is not so high and the drawing speed promote a significant cooling, the real fibre temperature is probably well approximated by the proposed correlation. Besides, the present correlation of the fibre temperature can be also utilized with the model of Buchak and Crowdy [36] that can handle fibres with general cross-sectional shapes.

5. Conclusions

In this work, we proposed a novel correlation to calculate the fibre temperature during the drawing process of PCFs. To this end, we assumed thermal equilibrium between the oven walls and the fibre. First, we utilized this correlation in the model of Chen et al. [32] to model the drawing process of an axis-symmetric capillary, and afterward, we used it with the model of Buchak et al. [34] to simulate the drawing process of holey fibres. The asymptotic models developed over the years by Stokes et al. [31], Chen et al. [32], and Buchak et al. [34] can be quickly solved numerically and does not necessitate significant computational resources compared to the full three-dimensional formulation. They can be solved on common workstations and the computational time is of the order of minutes which is significantly lower than the time required to carry out three-dimensional FEM simulations. The numerical predictions obtained in this work are overall in good accordance with experimental results despite the assumption that thermal equilibrium is only valid for fibre velocities close to the feeding speeds. A good concordance has been obtained both for the cases of SDRs and HDRs, both for annular capillaries and holey fibres, and both for pressurized and unpressurized fibres. The agreement between experiments and simulations could be even improved by selecting a different fibre temperature distribution, thereby furnishing experimenters with a powerful and fast predictive tool.