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] 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 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 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 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

Cumming and Howell [30], Stokes et al.[31], and Chen et al. [32] suggest the following scalings 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 while the dynamic and kinematic boundary conditions are transformed into and $ϵ$ and $Re$ are the slenderness parameter and the Reynolds number. They read

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$ 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 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 for $i=0,1,...N$. ${\overline{k}}_0^{\left(i\right)}$ assumes the form 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 Integrating now Eqs. (13a) and (13b) in the cross-sectional plane and applying the divergence theorem, one obtains after lengthy algebraic manipulations [31,32] where 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 Eq. (15) can be written in a Lagrangian reference frame where the problem is unsteady as 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).

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

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.

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).

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.