MicroStructure Element Method ( MSEM ) : viscous flow model for the virtual draw of microstructured optical fibers

We propose a new method to accurately model the structural evolution of a microstructured fiber (MOF) during its drawing process, given its initial preform structure and draw conditions. The method, applicable to a broad range of MOFs with high air-filling fraction and thin glass membranes, is an extension of the Discrete Element Method; it determines forces on the nodes in the microstructure to progressively update their position along the neck-down region, until the fiber reaches a final frozen state. The model is validated through simulation of 6 Hollow Core Photonic Band Gap Fibers (HC-PBGFs) and is shown to predict accurately the final fiber dimensions and cross-sectional distortions. The model is vastly more capable than other state of the art models and allows fast exploration of wide drawing parameter spaces, eliminating the need for expensive and time-consuming empirical parameter scans. ©2015 Optical Society of America OCIS codes: (000.4430) Numerical approximation and analysis; (060.2280) Fiber design and fabrication; (060.4005) Microstructured fibers; (060.5295) Photonic crystal fibers. References and links 1. J. C. Knight, T. A. Birks, P. S. Russell, and D. M. Atkin, “All-silica single-mode optical fiber with photonic crystal cladding,” Opt. Lett. 21(19), 1547–1549 (1996). 2. J. K. Ranka, R. S. Windeler, and A. J. Stentz, “Visible continuum generation in air-silica microstructure optical fibers with anomalous dispersion at 800 nm,” Opt. Lett. 25(1), 25–27 (2000). 3. F. Benabid, J. C. Knight, G. Antonopoulos, and P. S. J. Russell, “Stimulated raman scattering in hydrogen-filled hollow-core photonic crystal fiber,” Science 298(5592), 399–402 (2002). 4. P. S. J. Russell, P. Hölzer, W. Chang, A. Abdolvand, and J. Travers, “Hollow-core photonic crystal fibres for gas-based nonlinear optics,” Nat. Photonics 8(4), 278–286 (2014). 5. S.-C. Xue, R. Tanner, G. Barton, R. Lwin, M. Large, and L. Poladian, “Fabrication of microstructured optical fibers-part ii: Numerical modeling of steady-state draw process,” J. Lightwave Technol. 23(7), 2255–2266 (2005). 6. S. S. Chakravarthy and W. K. Chiu, “Boundary integral method for the evolution of slender viscous fibres containing holes in the cross-section,” J. Fluid Mech. 621, 155–182 (2009). 7. G. Luzi, P. Epple, M. Scharrer, K. Fujimoto, C. Rauh, and A. Delgado, “Numerical solution and experimental validation of the drawing process of six-hole optical fibers including the effects of inner pressure and surface tension,” J. Lightwave Technol. 30(9), 1306–1311 (2012). 8. E. N. Fokoua, M. N. Petrovich, N. K. Baddela, N. V. Wheeler, J. R. Hayes, F. Poletti, and D. J. Richardson, “Real-time prediction of structural and optical properties of hollow-core photonic bandgap fibers during fabrication,” Opt. Lett. 38(9), 1382–1384 (2013). 9. R. F. Cregan, B. J. Mangan, J. C. Knight, T. A. Birks, P. S. Russell, P. J. Roberts, and D. C. Allan, “Single-mode photonic band gap guidance of light in air,” Science 285(5433), 1537–1539 (1999). 10. C. M. Smith, N. Venkataraman, M. T. Gallagher, D. Müller, J. A. West, N. F. Borrelli, D. C. Allan, and K. W. Koch, “Low-loss hollow-core silica/air photonic bandgap fibre,” Nature 424(6949), 657–659 (2003). 11. P. J. Roberts, F. Couny, H. Sabert, B. J. Mangan, D. P. Williams, L. Farr, M. W. Mason, A. Tomlinson, T. A. Birks, J. C. Knight, and P. St J Russell, “Ultimate low loss of hollow-core photonic crystal fibres,” Opt. Express 13(1), 236–244 (2005). #225081 $15.00 USD Received 16 Oct 2014; revised 10 Dec 2014; accepted 12 Dec 2014; published 7 Jan 2015 © 2015 OSA 12 Jan 2015 | Vol. 23, No. 1 | DOI:10.1364/OE.23.000312 | OPTICS EXPRESS 312 12. F. Poletti, N. V. Wheeler, M. N. Petrovich, N. K. Baddela, E. Numkam Fokoua, J. R. Hayes, D. R. Gray, Z. Li, R. Slavik, and D. J. Richardson, “Towards high-capacity fibre-optic communications at the speed of light in vacuum,” Nat. Photonics 7(4), 279–284 (2013). 13. F. Poletti, M. N. Petrovich, and D. J. Richardson, “Hollow-core photonic bandgap fibers: Technology and applications,” Nanophotonics 2(5-6), 315–340 (2013). 14. D. G. Ouzounov, F. R. Ahmad, D. Müller, N. Venkataraman, M. T. Gallagher, M. G. Thomas, J. Silcox, K. W. Koch, and A. L. Gaeta, “Generation of megawatt optical solitons in hollow-core photonic band-gap fibers,” Science 301(5640), 1702–1704 (2003). 15. V. Sleiffer, Y. Jung, N. Baddela, J. Surof, M. Kuschnerov, V. Veljanovski, J. Hayes, N. Wheeler, E. Numkam Fokoua, J. Wooler, D. Gray, N. Wong, F. Parmigiani, S. Alam, M. Petrovich, F. Poletti, D. Richardson, and H. de Waardt, “High capacity mode-division multiplexed optical transmission in a novel 37-cell hollow-core photonic bandgap fiber,” J. Lightwave Technol. 32(4), 854–863 (2014). 16. A. M. Cubillas, J. M. Lazaro, O. M. Conde, M. N. Petrovich, and J. M. Lopez-Higuera, “Gas sensor based on photonic crystal fibres in the 2ν3 and ν2+ 2ν3 vibrational bands of methane,” Sensors (Basel) 9(8), 6261–6272 (2009). 17. S. Ghosh, A. R. Bhagwat, C. K. Renshaw, S. Goh, A. L. Gaeta, and B. J. Kirby, “Low-light-level optical interactions with rubidium vapor in a photonic band-gap fiber,” Phys. Rev. Lett. 97(2), 023603 (2006). 18. E. N. Fokoua, D. J. Richardson, and F. Poletti, “Impact of structural distortions on the performance of hollowcore photonic bandgap fibers,” Opt. Express 22(3), 2735–2744 (2014). 19. M. Matovich and J. Pearson, “Spinning a molten threadline. Steady-state isothermal viscous flows,” Ind. Eng. Chem. Fund. 8, 512–520 (1969). 20. S. R. Choudhury, “A computational method for generating the free-surface neck-down profile for glass flow in optical fiber drawing,” Numer. Heat Tr. Anal. Appl. 35, 1–24 (1999). 21. P. Gospodinov and A. Yarin, “Draw resonance of optical microcapillaries in non-isothermal drawing,” Int. J. Multiph. Flow 23(5), 967–976 (1997). 22. J. Ramos, “Drawing of annular liquid jets at low reynolds numbers,” Comput. Theor. Polym. Sci. 11(6), 429–443 (2001). 23. A. D. Fitt, K. Furusawa, T. M. Monro, C. P. Please, and D. J. Richardson, “The mathematical modelling of capillary drawing for holey fibre manufacture,” J. Eng. Math. 43(2/4), 201–227 (2002). 24. F. Geyling, “Basic fluid‐dynamic considerations in the drawing of optical fibers,” Bell Syst. Tech. J. 55(8), 1011–1056 (1976). 25. Y. Chen and T. A. Birks, “Predicting hole sizes after fibre drawing without knowing the viscosity,” Opt. Mater. Express 3(3), 346–356 (2013). 26. G. Luzi, P. Epple, C. Rauh, and A. Delgado, “Study of the effects of inner pressure and surface tension on the fibre drawing process with the aid of an analytical asymptotic fibre drawing model and the numerical solution of the full n.–st. Equations,” Arch. Appl. Mech. 83(11), 1607–1636 (2013). 27. R. Kostecki, H. Ebendorff-Heidepriem, S. C. Warren-Smith, and T. M. Monro, “Predicting the drawing conditions for microstructured optical fiber fabrication,” Opt. Mater. Express 4(1), 29–40 (2014). 28. Q. Sun, “Discrete modelling of two-dimensional liquid foams,” China Part. 1(5), 206–211 (2003). 29. Y. Maalej, M. El Ghezal, and I. Doghri, “Micromechanical approach for the behaviour of open cell foams,” Eur. J. Comp. Mech. 22, 198–208 (2013). 30. W. Warren and A. Kraynik, “The linear elastic properties of open-cell foams,” J. Appl. Mech. 55(2), 341–346 (1988). 31. Y. Kim and Y. Seol, “Numerical simulations of two-dimensional wet foam by the immersed boundary method,” Comput. Struc. 122, 259–269 (2013). 32. C. Pozrikidis, Boundary Integral and Singularity Methods for Linearized Viscous Flow (Cambridge University, 1992). 33. P. A. Cundall and O. D. Strack, “A discrete numerical model for granular assemblies,” Geotechnique 29(1), 47– 65 (1979). 34. G. Jasion, J. Shrimpton, Z. Li, and S. Yang, “On the bridging mechanism in vibration controlled dispensing of pharmaceutical powders from a micro hopper,” Powder Technol. 249, 24–37 (2013). 35. F. T. Trouton, “On the coefficient of viscous traction and its relation to that of viscosity,” Proc. R. Soc. Lond., A Contain. Pap. Math. Phys. Character 77(519), 426–440 (1906). 36. S. H. K. Lee and Y. Jaluria, “Simulation of the transport processes in the neck-down region of a furnace drawn optical fiber,” Int. J. Heat Mass Tran. 40(4), 843–856 (1997). 37. R. H. Doremus, “Viscosity of silica,” J. Appl. Phys. 92(12), 7619–7629 (2002). 38. Heraeus, “Heraeus quarzglas: Thermal properties, fused silica,” http://heraeusquarzglas.com/en/quarzglas/thermalproperties/Thermal_properties.aspx, Accessed 15/08, 2014. 39. Z. Wei, K.-M. Lee, S. W. Tchikanda, Z. Zhou, and S.-P. Hong, “Free surface flow in high speed fiber drawing with large-diameter glass preforms,” J. Heat Transfer 126(5), 713–722 (2004). 40. S.-C. Xue, L. Poladian, G. Barton, and M. Large, “Radiative heat transfer in preforms for microstructured optical fibres,” Int. J. Heat Mass Tran. 50(7-8), 1569–1576 (2007). 41. Y. M. Stokes, P. Buchak, D. G. Crowdy, and H. Ebendorff-Heidepriem, “Drawing of micro-structured fibres: Circular and non-circular tubes,” J. Fluid Mech. 755, 176–203 (2014). #225081 $15.00 USD Received 16 Oct 2014; revised 10 Dec 2014; accepted 12 Dec 2014; published 7 Jan 2015 © 2015 OSA 12 Jan 2015 | Vol. 23, No. 1 | DOI:10.1364/OE.23.000312 | OPTICS EXPRESS 313 42. Z. Yin and Y. Jaluria, “Neck down and thermally induced defects in high-speed optical fiber drawing,” J. Heat Transfer 122(2), 351–362 (2000). 43. J. Yang and Y. Jaluria, “Transport processes governing the drawing of a hollow optical fiber,” J. Heat Transfer 131(7), 072102 (2009).


Introduction
Since the first demonstration nearly 20 years ago [1] microstructured optical fibers (MOFs) have considerably diversified in topology, undergone an impressive improvement in performance and enabled novel applications, from the generation of ultra-bright supercontinuum light [2] to the demonstration of exciting light-matter interactions with unprecedented efficiency in compact gas-filled microstructured cells [3,4].
Optical performance improvements over this time have been driven by empirical trial and error in a fabrication process that is both costly and time consuming.As a result, even when the optimum structure to be targeted is known, development cycles of months, often years, are typically required to refine preform/cane structural parameters and identify drawing conditions to achieve it in practice.A numerical model that could predict the geometries of drawn fibers and highlight draw parameters necessary to achieve specific structures would be an invaluable tool for the further development of these fibers.However, the models which have emerged thus far are capable of simulating only very simple MOFs with a handful of holes [5][6][7][8].Here, we present an elegant physics-based model which can be used to accurately and efficiently model the fiber drawing process of MOFs with large air-filling fractions and thin membranes.Since the number of fibers that could be studied with this method is large -it comprises highly nonlinear MOFs as well as hollow core fibers guiding through antiresonance effects (e.g. with a full Kagome' lattice or with a simplified single ring cladding), or though photonic bandgap, we focus the presentation of the method on an example of a complex microstructured optical fiber, namely the hollow core photonic bandgap fiber (HC-PBGF).Comparison with experimentally fabricated HC-PBGFs demonstrates the accuracy with which the model can reproduce all the microstructured features of HC-PBGFs, both in terms of dimensions and of distortions.This model therefore provides a powerful tool to speed up fabrication and optimization of new microstructured fibers by providing the fabrication parameters necessary to achieve key microstructure features before a single fiber draw has been undertaken.Thus, it can reduce the uncertainty, expense and time necessary to produce high quality HC-PBGFs (or similarly structured MOFs).
HC-PBGFs are a type of optical fiber that uses an internal regular microstructure formed of hundreds of closely spaced holes to create a photonic bandgap that allows light guidance in an air core [9][10][11][12][13], Fig. 1.Since over 99% of the light is guided in air, HC-PBGFs are extremely resilient to non-linear effects, making them suited to high peak power delivery [14] and potentially to low latency high capacity data transmission [12,15].Besides, they can be used for gas for sensing [16], and efficient light-matter interactions [4,17].
The microstructure of a typical HC-PBGF is composed of slender glass membranes (~3μm x ~100nm) in a honeycomb arrangement with a central hollow core region (~10-30μm diameter), Fig. 1(c).This is normally fabricated in a 2 stage draw; an initial preform constructed of stacked capillary tubes (OD 2-4cm, Fig. 1(a)) is drawn through a vertical furnace to achieve transverse contraction and produce canes (OD 2-5mm, Fig. 1(b)).The cane is then placed in a thick jacket tube (OD 5-10mm) and drawn to optical fiber dimensions (OD ~125-200µm, Fig. 1(c)).In the second stage, as the assembled preform descends through the furnace, the temperature increases and the viscosity lowers dramatically; longitudinal elongation and diameter reduction are controlled through a capstan.During this process, surface tension drives a collapse of the holes in the structure which needs to be counterbalanced by an appropriate gas pressure applied by the fabricators to control the structure, Fig. 2. In this stage, fluid dynamics driven structural distortions can occur, resulting in fibers which are not a homothetically scaled version of the cane.Since the optical properties of the fiber depend on, and are extremely sensitive to, the geometry of each element of the microstructure and in particular those immediately surrounding the central core [18]; care must be taken to control it precisely to achieve a desired outcome.Fig. 2. Schematic of a second stage draw.Three regions are indicated; the preform which is the initial condition, the neck-down which is inside the furnace and is where all the geometry changes take place, and the fiber which is the final result.
The paper is organized as follows: in Section 2 we describe other model approaches from the literature.In Section 3 we describe the governing equations of our model and the numerical scheme.The parameters of our experimental and simulated fiber draws are reported in Section 4 and the results of these are discussed in Section 5.The conclusions are given in Section 6.

Existing models to study the drawing process of microstructured optical fibers
The modelling of microstructured optical fibers, of which HC-PBGFs are a specific kind, have roots in the initial models of solid fiber spinning [19,20], and have advanced through models of capillaries draws [21][22][23] to a small number of complete holey fiber models [5][6][7].
In our model we decouple the solution of the internal microstructure from the external fiber draw solution.In like fashion, we separate this background review into two sections: models for the external fiber draw process and methods for the internal microstructure solution.

Existing models for glass fiber drawing
The mathematics of the production of fibers using melt spinning has been considered for many decades.The focus of those studies has been to understand the neck-down shape that occurs as the large preform is drawn down to a slender fiber under heat and tension.An early mathematical approach can be found in [19], using a 1D approach from the Navier Stokes equations.The approach used in fiber spinning was applied to optical fibers, e.g [20,21,24], utilizing similar starting points but also including more advanced methods computationally and the effects of non-isothermal draws.
The distinguishing feature of most MOFs is the presence of holes that are continuous and uniform in the axial direction of the fiber.The theoretical study of capillary tube draws is therefore an appropriate starting point and has been performed in several variations both numerically and analytically [21-23, 25, 26] These models provide the neck-down shape of the internal and external diameter of a single tube, but cannot be easily extended to study multiple closely spaces holes.
A recent interesting approach, which employs similar ideas to those we use here, combines a capillary tube model with a nested microstructure [27]; here, a 3 hole microstructure is modelled within a capillary tube drawn using the Fitt model [23], with good agreement to experimental data.In our work we took a similar approach, but extended it for use with an HC-PBGF where hundreds of holes of arbitrary complexity are present.

Existing methods to model the microstructure evolution
The microstructure of an HC-PBGF is made up of glass membranes in a lattice like formation, Fig. 1(c).Methods have been developed to model similar structures, such as foam -both fluidic [28] and mechanical [29].The foam models aim to find homogenous bulk properties from analysis of a unit-cell [30], while these studies are not concerned with evolving specific elements of the structure the unit-cell analysis methods could be extended to resolve detailed deformation of the microstructure.The Immersed Boundary Method (IBM) was used by Kim and Seol [31] for a domain of arbitrary 2D foam.This method appears to be a potential option for modelling the HC-PBGF geometry, but to avoid the computational expense we opted for a simpler approach.
There have been at least 3 different approaches that look specifically at the evolution of the internal geometry of a microstructure fiber through a draw.Two of these use finite element (FE) methods, while the third uses the boundary integral method (BIM).In all cases however, only structures with a handful of holes were modeled [5][6][7].
Xue et al. [5], use the Finite Element method to draw a microstructure fiber with 4 holes.The results of Luzi et al. [7] demonstrate the use of FE for a 6 hole structure.The results of the study are in reasonable agreement with experiments, but it has yet to be demonstrated if these methods could be extended to a structure with hundreds of holes.The computational resources required were not specified in these studies, however, if mesh requirements scale linearly with strut number, the mesh size would be extremely high when applied to an HC-PBGF with ~100 × as many struts.
The Boundary Integral Method (BIM) is a numerical method that discretizes boundaries rather than volumes [32].This approach naturally lends itself to the slender elements that are less suited to FE meshes, and it is naturally designed for high viscosity flows.Chakravarty and Chie [6], use the BIM to model the evolution of holes in a fiber draw.In their simulations they tackle a system with 4 holes.Theoretically, the system could be extended to have many more holes, however, the 4 hole system took 100 hours to compute; for production use a faster method is clearly needed.Both FE and BIM methods resolve all the features of the surfaces, the methods operate at the microscopic level.So far, they have been applied to structures with a high volume to surface ratio; in contrast, the HC-PBGF structure, Fig. 1(c), has an extremely low volume to surface ratio and we can therefore expect FE and BIM methods to be too slow for practical uses.
The model we propose was inspired by the Discrete Element Method (DEM), but bares some relation to the beam methods used in the mechanical simulations of open cell foams [30].The DEM is a method that applies interaction forces to discrete bodies represented as Lagrangian particles and then evolves the system over time [33].DEM is primarily used to model granular flows and has demonstrated that a simple set of force expressions, when applied to large system of elements, is capable of recreating diverse and complex behaviors [34].A DEM methodology can be used in our rheology problem by considering the nodes as the Lagrangian particles.The pressure, viscosity and surface tension forces acting through the struts on those nodes are found, the node positions are updated and then the simulation is advanced in time.Compared with FE and BIM methods the accuracy with which surfaces are represented is cruder, however, the reduced computational resources enable the ability to represent genuine full-scale HC-PBGF geometries.
In this paper we develop and demonstrate the MicroStructure Element Method (MSEM), a numerical model to evolve a network of fluid membranes subjected to gas pressures, surface tension and viscous forces in a scenario that includes a non-isothermal furnace to represent the realistic scenario of MOF drawing.We apply this model to predict the final structure of HC-PBGFs; though the model would be suitable for other microstructure fibers composed of glass membranes, such as anti-resonant or suspended core fibers.Without loss of generality, we focus our attention and tailor the model to the study of HC-PBGFs with 7 rings of holes surrounding a core formed of 19 missing capillaries, and compare the simulated results with a set of fabricated fibers.

Numerical model
In this section we will describe the governing equations, structure and scheme of our numerical model.The numerical model presented here evolves a microstructured geometry through a first or second stage fiber draw and predicts the longitudinal evolution and final geometry of the resulting cane or fiber.The model has two parts addressing the microstructure and the jacket tube evolution.The microstructure is solved in 2D and stepped in the axial direction, z, from the initial preform through the furnace to the final structure.At each z position the microstructure is bound by the internal diameter of the jacket tube, h 1 .The neck-down shape of the jacket tube can be found using a model or measured data, in this work the neck-down shape is calculated using the expressions for a capillary draw developed by Fitt et al. [23].The Fitt model solves the neck-down geometry of the jacket glass independently from the solution of the microstructure geometry, Fig. 3.As the jacket tube contains significantly more glass than the microstructure geometry, this assumption is sensible and, as we shall see, sufficiently accurate.Through mass conservation the Fitt solution also finds the velocity and extensional change in the axial direction, which is passed to the microstructure model.

The microstructure model
The microstructure is represented in the model as a network of highly viscous fluid struts joined at nodes, Fig. 4(a).Figure 4(b) illustrates the forces acting on the struts around one node of the microstructure; surface tension is found at the nodes, and gas pressure and viscosity are found on the struts.These forces are used to evolve the microstructure through the draw, as explained in the following sections.connected by struts (black lines).(b) Forces: surface tension, F t , acts on the node directly and is a function of the angle turned in the fillet, ψ; the viscous, F v , and pressure, F p , forces act parallel and normal to the struts respectively.In this example the 3 left most nodes are adjacent to the core.The total force acting on the node is F n .

MSEM governing equations
Here we define the forces that act on the elements that make up the microstructure.HC-PBGFs are made from pure silica, which requires drawing these fibers at high temperatures, greater than 2000°C at the center of the hot zone in the furnace, where the viscosity of glass is μ ~10 5 Ns/m 2 .Using a typical strut length, l = 10μm, an estimate of strut strain rate, ε = 1μms −1 , and density, ρ = 2200 kg/m 3 , the resulting Reynolds number, Re = ρlε /μ~10 −14 , is extremely low indicating that the flow is entirely dominated by viscous forces rather than inertia.Thus a model that defines an internal force distribution on the nodes for a given boundary condition (the Fitt model) that is in equilibrium will be able to resolve the evolution of the internal structure.
We now consider the force balance on a single strut which we represent as a Lagrangian fluid volume, Fig. 5, with internal fluid forces dominated by viscous contributions.All mass in the microstructure model resides in the strut volumes, fluid is not permitted to flow from one volume to the next and the volume of each strut is conserved throughout the simulation.Two forces, per unit length, act on the strut volume, viscous stress and gas pressure.Fig. 5.A strut element represented as a Lagrangian fluid volume, the strut has length and width of l and w respectively and the direction parallel and normal to the strut is given by the unit vectors ŝ and n .P A and P B are the gas pressures on either side of the strut, and τ ss is the viscous stress, acting on the element.
The strut connects two nodes and the difference in velocity between these two nodes, u 1u 2 , generates a viscous stress, τ ss , in the direction parallel to the strut, indicated in Fig. 5 by the unit vector ˆ s .The viscous force per unit length, F v,s , acting on the strut is expressed by Eq. (1); subscript s indicates that the force acts on a strut, subscript n indicates that it acts on the node.

ˆ2 (
) where l and w are the length and width of the strut respectively, Fig. 5.The viscous force is a significant force throughout the draw, however, it does not drive geometry changes, it only acts to oppose changes to the structure and transmit changes across the structure.
A difference in gas pressure across the strut gives rise to the pressure force per unit length, F p,s , acting in the direction normal to the strut, n , given by Eq. ( 2).
( ) While active pressurization is used during fabrication, not every hole is controlled independently.Typically two different pressures are used during the draw: one to control the core and one for the cladding holes.In that case the pressure only acts directly on the struts around the core as there is no net contribution elsewhere in the cladding.This influences the position of the nodes on the core surround, and those changes are delivered to the rest of the cladding via the viscosity in the struts and surface tension at the nodes.
Surface tension acts on all the curved interfaces between the liquid glass and surrounding air.The general expression for the pressure caused by surface tension in 2D is given by Eq. (3).
where R c is the radius of curvature and γ is the surface tension parameter.In our 2D model we consider curvature in the axial direction to be negligible in comparison to the in-plane curvature, a reasonable assumption if one considers the ratio of length scales between the in plane features, i.e. strut length, versus the longitudinal features, i.e. furnace length, are of the order 10 −5 .
To calculate force due to surface tension a curved surface is assumed between two struts at their node, Fig. 4(b); in a first iteration we fit a circular arc in the junction.The surface tension force, F t,n,i , is the product of the area and surface tension pressure, Eq. (4)., , ˆ1 where the subscript i is the index for the junction, ψ is the angle turned in the junction, R c is the radius of the fillet and θ is the unit vector that is in the average direction of the two struts that make up the junction and it points away from the node, Fig. 4(b).The size of the fillet has no impact.The three forces: viscosity, pressure and surface tension; are illustrated in Fig. 4(b) with their directions relative to the struts indicated.
There is one surface tension force contribution for each strut-pair that is joined at that node, i.e. 2 or 3 contributions per node in the case of HC-PBGFs.The surface tension force acts at the node by definition.Net surface tension only exists on nodes whose struts are not uniformly distributed, i.e.where the struts intersect with the node at unequal angles.
To find the total force on each node the forces from the connecting struts are used, each node has 2 or 3 struts connected to it.The viscous, pressure and surface tension forces on a node are given by Eqs. ( 5)-( 7) respectively., , , , , , 1 2 , , , where the summation index identifies the struts joined to that node.The pressure force on a strut is halved as it is shared between its two nodes.The final step is to sum all the forces to calculate the total force, F n , acting on each node, Eq. ( 8), Fig.

The numerical scheme
The model is designed to produce a steady fiber geometry for any given set of steady draw conditions.The computational domain is discretized longitudinally and a non-isothermal scenario going from cold (outside the furnace) to hot (inside the furnace) and cold again is assumed (see Fig. 6).At any given step in z the forces are in equilibrium, therefore there should be no in-plane acceleration among the nodes in that z-slice; the numerical scheme finds the positions and velocities of the nodes that satisfy this condition.The scheme consists of 2 loops, the outer loop marches the solution from the start of the computational domain, z = 0, to the end of the computational domain, z = L, where the simulation stops.The inner loop finds the force equilibrium at each z step.

The outer loop
The outer loop simply advances z and enforces the external boundary of the microstructure to be the internal diameter from the Fitt model, h 1 (z) (see below).Before the inner loop starts the system suggests the final position and velocity of the nodes in this z slice to speed up convergence.Velocity is estimated using: u 1 = u 0 , the position is estimated as , where superscripts 1 and 0 refer to current and previous z steps.

The inner loop
For any given z slice the inner loop performs the following tasks: 2. The positions and velocities are advanced in the direction of the total force using Eqs.( 9) and (10), where a is the under relaxation factor which controls how quickly changes are introduced, and Δ Δ / dz dt t z = from the Fitt model.u u aF = + (9) 3. The strut lengths are dictated by the distance between the nodes they connect, the width is updated at each iteration and conserves the mass of the strut (also accounting for longitudinal elongations).
4. After the positions have been updated the convergence criterion is tested, Eq. ( 11).The inner loop then iterates until such criterion is satisfied.We find C = 1 × 10 −4 to be a good trade-off between accuracy and speed.
( ) 5. Once convergence has been achieved the final positions and velocities at this stage are stored, and the solution advances in z with the outer loop.For each solved z slice the stored node positions and velocities are also used for the prediction of the next z slice.In this manner the positions of the nodes advance in z; the velocities, however, are used to suggest the likely solution of the next z slice -this speeds up convergence but is not necessary and so we consider the velocities (in this implementation) to be independent from one slice to the next.

Numerical parameters
There are several numerical parameters that need to be defined for the effective implementation of this model.To appropriately resolve the changes during the draw ∆z must not be too large, we found all values of ∆z ≤ 10 −3 to give final geometries that differ from each other by no more than 0.01%, in this investigation.
The under relaxation factor, a, in Eq. ( 9) must be small enough to avoid divergence but large enough to allow convergence in a reasonable amount of time.The computation time of the scheme scales linearly with the number of nodes, there are ~700 nodes in this set of simulations and simulations take ~10 minutes.

Initial conditions
At z = 0, where the fiber is still 'cold', all node velocities are zero, and the positions and widths are prescribed to match the preform geometry.

Boundary conditions
During the simulation the nodes at the external limit of the cladding are bound to the inner diameter of the jacket tube, h 1 , defined by the Fitt solution, their velocity is fixed at zero, Eq. ( 12) and ( 13).The contraction rate of h 1 is taken into account by the first term on the right hand side of Eq. ( 10) in this manner the boundary contraction rate is removed from the active part of the scheme while the contraction rate contribution is fully incorporated into the solution.
where subscript ext indicates the indices of nodes at the external boundary of the microstructure, 0 ˆext x is the unit vector of the previous iteration external nodes and B is a factor relating to the protrusion of packing rod glass.The structure is centered at the origin.All other nodes are free to move as the scheme dictates.When a physical cane is made the gaps around the final ring of tubes are filled with packing rods.The packing rod glass influences the shape of the microstructure boundary which has been included here in the factor B, Eq. ( 14), which is a function of the expansion ratio, e, Eq. ( 15).ext,cane 1 2 10 10 20 where R ext,cane is the radial distance of the external nodes, from the centre, when originally defined at the start of the simulation.

The Fitt model for the bounding jacket glass
To apply the microstructure element model to the task of simulating an HC-PBGF draw we use a capillary draw to define the external boundary of the microstructure through the neckdown region.The Fitt model [23], is a 2D axisymmetric model, derived from the Navier-Stokes equations, to solve the shape of a glass capillary draw.Equations ( 16)-( 18) are from the Fitt model and are used to find the neck-down shape, the reader is referred to the original work for a detailed derivation; Eq. ( 19) is the energy equation based on the original Fitt model, however, we have adjusted the original to avoid a units discrepancy in the radiation term.We assume, following Fitt, that only viscosity, μ, is a function of temperature and that liquid silica is an incompressible Newtonian fluid.We also assume that the initial geometry of the tube is known along with feed speed and the temperature profile of the furnace, T a .Equations ( 16)-( 19) define the inner, h 1 , and outer, h 2 , diameter, axial velocity, w, and temperature, T, as functions of z.

t z p h h h h h h h h w h h
( ) ( )

t z p h h h h h h h h w h h
where subscripts of t and z indicate derivatives in time and the axial direction, g is acceleration due to gravity, γ and ρ are surface tension and density of the glass, p 0 is the gas pressure inside the tube relative to ambient pressure.In Eq. ( 19) T is temperature, k is thermal conductivity, σ is the Steffan Boltzman constant, α is emissivity and N is a heat transfer coefficient.
The left hand term in Eq. ( 16) is the inertial term, on the right hand side the terms relate to the viscous force and the surface tension force in the axial direction.Equations ( 17) and (18) apply the pressure which competes directly with surface tension, the rate of which is governed by the viscosity.Equation ( 19) solves the temperature of the glass, the terms on the left hand side are thermal inertia, conduction and radiation, the term on the right hand side relates to cooling.We specify the entire furnace temperature profile, which we assume to have a Gaussian longitudinal profile in the furnace region and a step change to ambient temperature at the end of the furnace, the domain continues until the fiber has cooled and no further changes occur.The glass temperature at z = 0 we assume to be at the same temperature as the furnace at z = 0 as the preform is moving very slowly here and has been in the vicinity of the top of the furnace for some time.The thermal conduction term is small and we choose to neglect it.
We are interested in the steady state structure so assume all temporal derivatives to be zero.The resulting expressions are solved using a simple Euler scheme marching forward axially through the neck-down region from z = 0 to z = L.The solution to the Fitt model equations is found using the shooting method targeting the final fiber OD.The pressure applied in the tube can also be included in the shooting method to target a specific expansion ratio, Eq. ( 15).
Solving Eqs. ( 16)-( 19) results in a complete neck-down shape for a given set of draw parameters, an example is given in Fig. 6.Comparing the mass flow rates at the beginning and end of the domain gives ( ) Fig. 6.A neck down profile generated by the Fitt model (bottom) and temperature profile of furnace and preform / fiber (top).The draw diameter and draw speed are relative to the initial diameter and feed speed.The 'node transect' shows the positions of a selection of nodes from the microstructure as they are evolved through the draw, (data from simulated fiber 1 -see Table 3).

Fiber draws
In this section we describe the fiber draw experiment that we have conducted to validate the model.We drew 6 different HC-PBGFs from the same cane and independently simulated them with the MSEM.Note that although the model is also capable of simulating cane draws from first stage preforms, since no active pressure is applied and consequently the observed deformations are typically much less pronounced in this case, the results are visually less interesting and we do not include them here.Here the model is validated by simulating the draw of cane into fibers.

Experimental draw
A 19 cell 7 ring cane, 3mm in diameter, was placed in a jacket tube with inner and outer diameter of 3.5 and 10.0 mm.The cane can be seen in Fig. 7(a).In total 6 fibers were drawn, with 2 nominal expansion ratios e (defined as the microstructure diameter to outer diameter ratio, akin to Eq. ( 15)), each with 3 different core pressure values but constant cladding pressure.Fibers A-C had a nominal expansion ratio of e ~50% and fibers D-F had an expansion of e ~60%.The pressure of importance in this study is the difference between the cladding and the core pressure, ∆P core = P core -P cladding .The pressure of Fiber A was chosen to produce a typical sized core, this became the baseline pressure.The pressures of the following 2 fibers were 4 kPa and 8 kPa above such baseline.All other parameters were kept constant.The draw conditions can be seen in Table 1.Such a wide range of pressures was chosen because it generated fibers with distortions that exceed those normally considered acceptable, which allowed us to test our model even beyond what would be its normal operating regime.

Virtual draw (MSEM simulations)
To describe the MSEM simulations that virtually reproduce the real draws discussed above we must first describe in more detail our initial conditions, material parameters and furnace assumptions.
The initial condition is the cane geometry used in the experiment, combined with the jacket glass.We assume there is no void between the cane and jacket glass and instead extend the cladding around the cane so that the total glass volume is matched.Figure 7 compares the real cane with its MSEM representation.

Furnace profile
The simulations are non-isothermal, and choice of cold enough initial and final steps provides an explicit beginning and end: before changes occur and after all changes have ceased.The longitudinal temperature profile of the fiber furnace at the Optoelectronics Research Centre was measured with a thermocouple at several peak temperatures and was found to be approximately Gaussian in the axial direction, z, along the central axis.The computational domain starts at the top of the furnace, proceeds through the peak temperature region (the 'hot zone') to the furnace exit, and includes a sufficiently long region of space below it to allow the fiber to cool.The fiber is considered cool when the outer diameter stops changing, dh 2 /dt ≈0, Fig. 6 shows the fiber speed asymptote as the diameter becomes steady.
For these simulations the length of our computational domain is L = 0.1 m, the furnace region is 0.08m and starts at z = 0, the peak temperature is T peak = 1822 °C at z = 0.04 m with a standard deviation of σ = 0.06 m and an ambient temperature of 20°C, Fig. 6.The fiber temperature is found by solving Eqs. ( 21)-( 24) from the Fitt model.The peak temperature was chosen such that simulated fiber tension,  19) was found to have a significant effect on the temperature profile and neck down shape [36], however, if tension values were matched, the microstructure geometry was found to be insensitive to changes of N.

Material parameters
The HC-PBGFs developed at Southampton are made from pure silica, Heraeus F300 glass.
For the temperature dependence of its viscosity we used the fitted curve reported in [37], Eq. ( where T has units of °C, μ is in Ns/m 2 and R is the universal gas constant.The remainder of the materials parameters are given in Table 2.

Virtual draw parameters
The draw parameters of the virtual draws are given in Table 3. Experimental draws are labelled A-F, while virtual draws (MSEM simulations) are labelled 1-6.3).Left: the axial temperature profile of fiber and furnace.Center: complete neck down shape across the whole computational domain with initial and final geometries shown above and below respectively.Right: the microstructure at specified positions in the neck down, with fiber OD and cladding expansion ratio at that position, and relative force vectors indicated.Microstructure changes are concentrated between z = 10 and 50%.
The ∆P core used for fiber 1 was chosen to produce the same core size as that of experimental fiber A; fibers 2-6 use this pressure baseline and then step 4 kPa and 8 kPa as in the experimental draw.All the pressure values in the simulation are offset by 1.8 kPa relative to their equivalent experiment, the origin of this offset is discussed later.

Results
As shown below, one immediate benefit of the MSEM results is that they allow us to study longitudinal structural evolutions at any point inside the furnace and understand where a given force dominates, or identify the length scale for changes to occur.
Figure 8 shows the neck down shape defined by the Fitt model and shows how the microstructure evolves as it is drawn down.The extent of the computational domain is chosen to contain all the fiber changes, however, the majority of the changes occur in the hot zone, between z = 10 and 50%.The first change is in extension; the OD of the fiber reduces rapidly while the ID remains constant -the result is a thinning of the microstructure struts, Fig. 8(a,b).The struts continue to thin almost right through to the end.The pressure force dominates at the start and causes considerable deformation around the core, Fig. 8(b).As the geometry reduces in size the pressure force becomes less significant; by 50%, Fig. 8(c), the surface tension is influencing the structural changes more than the pressure and the core begins to enlarge.As the glass cools the viscosity increases and counters further changes by the surface tension, Fig. 8(d).By 50%, Fig. 8(c), the shape of the microstructure is close to its final geometry, Fig. 8(d); very little change is distinguishable between this stage and the final fiber save for a reduction in OD.
It was found that the length of the computational domain required is based on the rate of temperature change of the fiber, Fig. 6.If the fiber heated and cooled more rapidly the viscosity would increase more rapidly and there would be less axial length within which changes could occur.However, additional simulations (not presented here) indicate that increasing the cooling rate while continuing to target the same final OD, expansion ratio, and tension, results in an almost indistinguishable microstructure.The analytical works of Chen and Birks [25], and Stokes et al. [41] found that the final geometry of a drawn hole or hollow fiber can be predicted by knowing the tension, and that the neck down shape is not required.It is surprising, however, to find this the case for the core size with separate pressure control.This brief test, however, is by no means conclusive.
The simulated fibers are compared with the experimental fibers qualitatively, examining the structural deformations, and then quantitatively by comparing the core size for each pressure.Figure 9 compares optical micrographs of experimental fibers A, B & C with the corresponding simulated fibers 1, 2 & 3. First consider the core size; both experiments and simulations follow the trend of expanding core size with increasing core pressure.The general cladding structure is reasonably uniform in shape and size when the core is small, but as the core is expanded the cells in the surrounding structure begin to compress.The compression is felt strongest in the ring nearest to the core while the cells near the jacket glass show the least (but still notable amounts of) deformation.Finally, examine the first ring nearest to the core, the corner cells are larger and have a different shape relative to the side holes in between them.As the core expands the corner holes stretch until, in the final case, the membranes are so thin they are unobservable with an optical microscope.All of the features seen in the experimental micrographs are recreated in the simulation with apparently equal magnitude.
Figure 10 gives a quantitative comparison between all the experiments and simulations using the normalized core size, R core / R clad , the values of which are reported in Tables 1 and 3.The core size is plotted against pressure offset from fiber A or 1 for experimental and simulated results respectively, as indicated in Tables 1 and 3. Figure 10 shows an impressive match between the experiment and the simulation.The difference between core ratio of the experiments and the simulations does not exceed 4%.The range of core sizes considered encompasses the desirable core size bracket of these types of fibers and demonstrates that the numerical model is well bounded.While the agreement with fabricated fibers was overall very encouraging, there were nonetheless some discrepancies in the drawing parameters which are worth commenting on.
First of all, there was a constant 1.8 kPa pressure difference between the experimental and simulated fibers.While the source of this difference is still to be exactly pinpointed, we can suggest possible explanations.Since the experiments required a larger pressure difference it is possible that some leak exists in the gas flow system.After the pressurized gas has been measured, it flows down narrow pipes and glass channels, and it is well known that flow through narrow pipes will add a pressure drop to any gas flow.Another potential source of the discrepancy is that the model assumes a uniform radial temperature profile, work by Xue et al. [40] is just one of a number of studies considering non-uniform radial temperature profiles.If the glass in the vicinity of the core was cooler then the viscosity would be greater and more pressure would be required to cause core expansion.
Moreover, another discrepancy we found was that in order to ensure that the simulated draw tension matched the experiment, we had to impose a cooler furnace in our simulations than in the experiments.Since the draw tension is indicative of viscosity of the fiber, we believe this is a more appropriate parameter to match than the outside temperature.Fitt et al. [23] admit the thermal model and boundary conditions may require more detail and the heat transfer coefficient is difficult to estimate [42] and could change axially [39] as well as radially [43].These are all areas for future investigations.

Conclusions
A microstructure element model (MSEM) that predicts the structural evolution of a hollow core photonic band-gap fiber during its drawing process has been developed, demonstrated and validated.The model combines the capillary draw expressions of Fitt et al. with a bespoke Lagrangian fluid element model that solves viscous fluid forces across the microstructure.We tested the model by virtually drawing a range of fibers from the same initial condition and compared the results with an experiment of the same.A constant pressure offset was applied to all the simulations, in such a way to match the core size of the first simulation with the first experiment.With the offset in place all the other results match remarkably well.The model, while simple in conception, demonstrates a capability that is far beyond the previous state of the art.While the model presented here has been applied to a specific HC-PBGF, it must be stressed, that it can be easily extended to any other HC-PBGF design and indeed to many other MOFs that are composed of thin glass struts.Suspended core and anti-resonant fibers are two examples that the MSEM could easily be adapted to solve.Likewise, non-Newtonian polymers could also be simulated by replacing the viscous force, Eq. ( 1), with a viscoelastic force expression.With a typical simulation time of ~10 minutes for a full virtual draw, fast exploration of vast drawing parameter spaces is possible.This can eliminate the need for expensive and time-consuming empirical parameter scans and has the potential to considerably speed-up future development of high performance MOFs of different types.

Fig. 3 .
Fig. 3.The relationship between the microstructure solution and the Fitt solution of the Jacket glass.The external boundary of the microstructure solution is bound by the internal diameter of the jacket, h 1 , as it marches in z.The jacket tube is solved using the Fitt model.

Fig. 4 .
Fig. 4. Geometry representation and force resolution around a node.(a) The microstructure is represented by a triangular lattice of cells.This generates a grid of nodes (red points)connected by struts (black lines).(b) Forces: surface tension, F t , acts on the node directly and is a function of the angle turned in the fillet, ψ; the viscous, F v , and pressure, F p , forces act parallel and normal to the struts respectively.In this example the 3 left most nodes are adjacent to the core.The total force acting on the node is F n .

Fig. 7 .
Fig. 7.The cane as used in the experiments, (a), and simulation, (b), without the jacket glass.

Fig. 8 .
Fig. 8. Evolution of the microstructure through simulated draw 1 (Table3).Left: the axial temperature profile of fiber and furnace.Center: complete neck down shape across the whole computational domain with initial and final geometries shown above and below respectively.Right: the microstructure at specified positions in the neck down, with fiber OD and cladding expansion ratio at that position, and relative force vectors indicated.Microstructure changes are concentrated between z = 10 and 50%.

Fig. 10 .
Fig. 10.Core size as a result of different core pressures, comparing the experimental results with the simulated results of all fibers.