Unraveling radial dependency effects in fiber thermal drawing

Fiber-based devices with advanced functionalities are emerging as promising solutions for various applications in flexible electronics and bioengineering. Multimaterial thermal drawing, in particular, has attracted strong interest for its ability to generate fibers with complex architectures. Thus far, however, the understanding of its fluid dynamics has only been applied to single material preforms for which higher order effects, such as the radial dependency of the axial velocity, could be neglected. With complex multimaterial preforms, such effects must be taken into account, as they can affect the architecture and the functional properties of the resulting fiber device. Here, we propose a versatile model of the thermal drawing of fibers, which takes into account a radially varying axial velocity. Unlike the commonly used cross section averaged approach, our model is capable of predicting radial variations of functional properties caused by the deformation during drawing. This is demonstrated for two effects observed, namely, by unraveling the deformation of initially straight, transversal lines in the preform and the dependence on the draw ratio and radial position of the in-fiber electrical conductivity of polymer nanocomposites, an important class of materials for emerging fiber devices. This work sets a thus far missing theoretical and practical understanding of multimaterial fiber processing to better engineer advanced fibers and textiles for sensing, health care, robotics, or bioengineering applications.

Fiber-based devices with advanced functionalities are emerging as promising solutions for various applications in flexible electronics and bioengineering. Multi-material thermal drawing in particular has attracted a strong interest for its ability to generate fibers with complex architectures. Thus far however, the understanding of its fluid dynamics has only been applied to single material preforms for which higher order effects, such as the radial dependency of the axial velocity, could be neglected. With complex multi-material preforms, such effects must be taken into account, as they can affect the architecture and the functional properties of the resulting fiber device. Here, we propose a versatile model of the thermal drawing of fibers which takes into account a radially varying axial velocity. Unlike the commonly used cross-section averaged approach, our model is capable of predicting radial variations of functional properties caused by the deformation during drawing. This is demonstrated for two effects observed, namely by unraveling the deformation of initially straight, transversal lines in the preform and the dependence on the draw ratio and radial position of the in-fiber electrical conductivity of polymer nano-composites, an important class of materials for emerging fiber devices. This work sets a thus far missing theoretical and practical understanding of multi-material fiber processing to better engineer advanced fibers and textiles for sensing, health care, robotics or bioengineering applications.
Thermal drawing is at the heart of the fabrication of telecommunication optical fibers. As illustrated in Fig. 1(a,b), it consists in heating a macroscale preform, typically made out of glass or thermoplastics, to its softening temperature and pulling it to create a much thinner and longer fiber that maintains the initial cross-sectional architecture. From simple glass-based step index structures, the design of optical fibers was expanded to single-mode waveguiding, photonic crystal and Bragg mirror fibers [1][2][3][4][5][6][7] . This technique has also been used to make micro-structured multi-material fibers that integrate optical, electronic and optoelectronic materials [7][8][9][10][11] . These fibers exhibit advanced functionalities 9,12 , making them promising building blocks for applications in optics but also soft electronics [13][14][15][16] , optoelectronics 10,11,13,17,18 , sensing 16,19 or bioengineering 19-21 . While the introduction of different materials in the drawing process has triggered much interest to realize innovative fiber-based devices and smart textiles, subtle effects of the fluid dynamics of multi-material co-drawing remain to be investigated to better understand and exploit this approach. In particular, modeling and analytical analysis have mostly relied on fluid dynamics descriptions that assume a velocity in the drawing direction that is independent of the radial position 19,22,23 . There exist some numerical studies that are not based on such an assumption, but only in the context of single material, micro-structured optical fibers 23,24 . The impact of the non-uniform velocity field in the cross-section in the context of multi-material drawing is however essential, albeit yet unexplored. A radially varying velocity field and hence deformation rate during drawing can influence the targeted architecture and resulting properties of the functional a) Alexis G. Page and Mathias Bechert contributed equally to this work. b) Electronic mail: fabien.sorin@epfl.ch FIG. 1. (a) Visualization of the thermal drawing of a fiber in a furnace with three heating zones. We combined a picture of an interrupted draw with a schematic of the furnace. (b) Sketch of the axisymmetric model. The ambient temperature T air is assumed to have a Gaussian shape along the z-axis. materials at the fiber level. A striking manifestation of this effect is schematically shown in Fig. 1(a), where initially rectangular domains of polymer nano-composites deform significantly during drawing. The pattern observed, which we experimentally show and model below, results from a relative displacement of material at different radial positions. Taking this effect into account is primordial to realize advanced fiberbased devices with controlled properties.
In this work, we report an in-depth fluid dynamics analysis of the thermal drawing process that takes into account the radial dependency of the axial velocity. We propose a model that includes higher order terms of the velocity field, which enables us to elucidate the observation of line deformation shown in Fig. 1(a). From this velocity field, we can extract a more accurate distribution of the rate of deformation during drawing, which is crucial for many material properties at the fiber level. In particular, we apply the results of our model to precisely account for the dependence of electrical conductivity of polymer nano-composites on the radial position and the draw ratio, i.e. the ratio of preform to fiber diameters. Conductive polymer composites form an important class of materials widely used in fiber-based devices for optoelectronics 10,17,18 , bioengineering 20,21 , and electromechanical sensing 15,25 . They however have a well-known, yet so far unexplained dependence of their electrical conductivity on the fiber drawing conditions 14 . Our model advances the understanding of multi-material fiber drawing and opens opportunities for the design of advanced fiber-based devices and smart textiles with optimum control over material microstructure.
We extend the reduced one-dimensional model for elongating fibers widely used in literature [26][27][28][29] to cover thermal effects and to access information about the radial distribution of the variables. We consider only the cladding material in the model, as the functional materials represent only a small fraction of the fiber cross-sectional area, and assume an axisymmetric, round fiber. A detailed derivation of the equations is given in the supplementary material. As shown in Fig. 1(b), the fiber can be described by its radius R(z) along the main axis together with the flow velocity field (u, v), with u(r, z) and v(r, z) denoting the radial and axial velocity components, respectively, as well as temperature T (r, z) and shear viscosity η(r, z). All variables are made dimensionless using the initial feeding speed v 0 , the drawing length L, maximum air temperature T max and corresponding minimum viscosity η min as reference values. Assuming a slender fiber 29,30 , i.e. α 2 = (R 0 /L) 2 1, enables us to expand the variables in α 2 , e.g., with the superscript indicating the order of expansion. Using an effective drawing length of L = 40 cm and an initial radius of R 0 = 1.25 cm, α 2 = O(10 −3 ). Following a previous study 31 , we assume a parabolic temperature field in r, which enables us to calculate the cross-section averaged governing equations. The cross-section average is denoted by · , and the Biot number, which compares the heat exchange between the material and the surrounding air to the heat conduction within the fiber, is defined by Bi = R 0 h/k, h being the heat transfer coefficient and k the material heat conductivity. Note that Eq. (2) implicitly assumes that the preform is primarily heated by the surrounding air, even though weak radiative heating can be accounted for by using this approach as well 32 . In some thermal drawing processes radiative heating can lead to a material temperature higher than the air temperature 33 , which cannot be modeled by this approach. Similarly, some functional materials may change the mechanical and heat properties of the fiber significantly, despite their low proportion inside the fiber. For the description of the surrounding air temperature, we assume a centered Gaussian profile, with dimensionless amplitude Λ and width ∆.
Utilizing the expansion in α 2 , the steady state leading order equations are obtained by cross-section averaging the continuity, momentum and heat equations, which yields with radially independent axial velocity v (0) = v (0) (z). The Stanton number St compares the heat transfer between the fiber and the surrounding air to the heat advected by the polymer flow and is given by St = 2 Bi/(α Pe), with the Peclet number Pe = ρ c p v 0 R 0 /k, c p denoting the heat capacity of the material. The energy equation (4c) is coupled to the momentum equation (4b) using Arrhenius' law, with dimensionless activation energy µ. The streamwise boundary conditions are given by the process setup as with the draw ratio Dr defined by the ratio of preform to fiber diameters. h, Λ and ∆ are determined using the experimentally determined preform-to-fiber radius profile during drawing. A detailed description of the parameter determination procedure is given in the supplementary material. We solve Eqs (4-6), which will also be referred to as the cross-section averaged model in the following, by numerical continuation using AUTO97 34 to obtain R, v (0) , u (0) = −r ∂ z v (0) /2, T (0) and η (0) . Using these solutions, we then go beyond the usual analysis by calculating the radially varying T (0) and η (0) , as well as investigating the α 2 -order of the expanded governing equations, which yield Equation (7) can be solved by numerical integration to obtain v (1) and u (1) = r 0 dr (r ∂ z v (1) )/r, both depending on r and z. To demonstrate the importance of this last step, we calculate the evolution of an initially straight line perpendicular to the drawing direction during the drawing process.  , an initially straight line remains straight during the entire process, which is not surprising as the dominating axial velocity does not vary along the radius at leading order. However, if the higher order terms u (1) and v (1) are included in the calculation (red lines), the lines clearly deform during the drawing, finally looking like an inverted 'V'. In order to visualize this effect experimentally, we fabricate (see supplementary material) a polysulfone (PSu) preform of rectangular cross-section (2.5 cm × 1.4 cm) within which 1 mm-wide sheets of carbon-black-loaded polycarbonate (PC) are equidistantly placed, perpendicular to the drawing direction. The feeding speed is set to 2 mm/s and the draw ratio is Dr = 8.31. The furnace has three heating zones (see Fig. 1(a)), which are set to 200/300/120 • C. The drawing of this preform is interrupted so that the sheets visualize the deformation of initially straight, horizontal material lines during drawing, as shown in Fig. 2(b).
The calculated (red) lines coincide very well with the ob-served shape of the sheets in the PSu fiber. Note that this comparison is nevertheless of qualitative nature, as the fiber used in the experiment has a rectangular cross-section, while the model assumes a circular fiber. A PC fiber is found to yield less quantitative, but still qualitative coincidence with the model prediction, which is probably due to the simplification of heating mechanisms as discussed above. For both materials, however, the final shape of the deformed lines is reproduced with high accuracy. All parameter values are given in the caption of Fig. 2. Physically, the non-uniform propagation of fluid elements lying on an initially straight line is the result of an interplay between a radially varying viscosity (see Fig. 2(a) and supplementary material) and the deformation of the fiber due to drawing. It is worth mentioning that the significant influence of the higher order velocity terms is a cumulative effect. At one particular position in z, the change in velocity introduced by the higher order is at most only a few percent. This is expected, as the contribution is scaled by α 2 according to Eq. (1), with α 2 being small at the basis of the expansion. By following a material point during the drawing, however, we integrate along the fiber length, and by that the order of magnitude increases from α 2 = (R 0 /L) 2 to R 2 0 /L. For this reason, the widely used cross-section averaged model, Eqs. (4)(5)(6), is well-suited if one is interested exclusively in the fiber radius and the averaged velocity field, as it is the case for single-material fibers, but insufficient if radially distributed properties come into play, like for multi-material fibers.
To further highlight the importance of this understanding, we now investigate the influence of thermal drawing on the electrical conductivity of nano-composites. As conductive fillers, we use carbon black (CB) and carbon nanotubes (CNT). In particular, we fabricate two types of preform: a PSu cladding with a PC-CNT composite sheet and a PC cladding with a PC-CB composite sheet. For drawing the preform with PSu cladding, we use the same settings as above. The preform with PC cladding has a rectangular cross-section (2.5 cm × 2.3 cm), the feeding speed is set to 2 mm/s, and the heating zones of the furnace are set to 120/270/120 • C.
As visible in Figs. 3(a,b), the nano-composite sheet is in contact with several metallic ribbons (2 mm × 1 mm crosssection) made of a bismuth-tin (tin-zinc) alloy for the PC (PSu) cladding. Even though the ribbons melt during consolidation and drawing, the high viscosities of the cladding and sheet materials ensure that the cross-sectional architecture is preserved. By polishing away the cladding of the drawn fiber, we can individually contact the electrodes with a multimeter to measure the conductivity perpendicular to the drawing direction (see Fig. 3(a)).
In order to access theoretically the final conductivity, we use a simple percolation law to link the effective filler concentration p e to the conductivity σ , with initial conductivity σ in of the undeformed material having an initial filler concentration p in , critical concentration p c , and critical exponent c. To account for the destruction of conductive pathways due to deformation during drawing we propose an empirical kinetic equation inspired by a previous work 35 , which links the effective filler concentration to the second invariant of the deformation tensor √ II D (see supplementary material), a denotes a destruction rate and analogue to the previous work 35 , we set m = 1.5. Solving Eq. (9) and substituting in (8) then yields with γ(r f ) denoting the trajectory of a material element during drawing with final radial position r f . The product a c and σ in are determined by fitting the experimental data. Equation (10) yields the radial distribution of conductivity in the drawn fiber, which can then be used to calculate the averaged conductivities measured between two electrodes (see Fig. 3(a)). Fig. 3(c) shows the perpendicular conductivity of the PC-CNT composite in a PSu fiber, measured between two electrodes placed at r/R = ±0.6 for a large range of draw ratios, and compares it to the theoretically determined averaged conductivity (red line). The increasing deformation with increasing draw ratio leads to a change in conductivity over several orders of magnitude, which is well reproduced by our model. Also plotted is the theoretical prediction based on the crosssection averaged model (black dashed line). Both models lead to practically identical results, which is not surprising as we average the conductivity over a large part of the fiber radius. In order to access experimentally the conductivity as a function of the radial position, we prepare a PC/PC-CB preform with six electrodes placed at equal distances along the fiber width (see Fig. 3(b)). Fig. 4 shows the measured conductivities for two different draw ratios. The points are placed in the middle of the domain of average, which has a size of 2.5 mm, i.e. 0.25 R, corresponding to the distance between the electrodes. These measurements are compared to the theoretically obtained moving-averaged distribution (red line). Despite its simple approach, our model describes the variation of conductivity along the radial position very well. Note that analogue to the PSu/PC-CNT fiber, the values of σ in and a c are uniquely set, independent of the draw ratio. Also shown in Fig. 4 is the averaged conductivity distribution resulting from the cross-section averaged model (black dashed line). Clearly, this model is incapable of predicting such a large radial variation of the conductivity, emphasizing the need of our extended model. The model presented here is able to predict the radial distribution of material properties in multimaterial fibers fabricated by weakly radiative thermal drawing, while it maintains the simplicity and flexibility of the cross-section averaged onedimensional model for single-material fibers. Even though assumptions like using a circular cross-section or employing the simple empirical kinetic equation (9) appear to be rather strong, we are able to reach qualitative and quantitative coincidence with our experimental findings. The flexibility of the model enables comprehensive parametric studies to find routes to innovative experimental designs. Our approach improves the understanding of multi-material fiber drawing and opens opportunities for the design of advanced fiber-based devices with optimum control over material microstructure and tailored properties.

SUPPLEMENTARY MATERIAL
The supplemental materials provide a detailed derivation of the thermal drawing model, a description of the materials and experimental methods used in this work, all parameter values used in the modeling, together with a description of how they are obtained, and further information on the evolution of material temperature, viscosity and velocity fields during drawing.

Governing equations
We follow the derivation as presented by Bechert and Scheid [1] for fibers and extend the model equations to cover thermal effects, analogue to Scheid et al. [2]. Moreover, we analyze higher order terms of the velocity field.
We model the fiber material as an incompressible, Newtonian fluid and assume negligible effects of inertia, gravity and surface tension (Re < 10 −6 , Re/F r < 1, Ca > 10 6 ). We further neglect all time derivatives as we are interested in steady solutions only. The corresponding governing equations are the continuity, momentum and heat equations, which can be written as, respectively, with ∇ = e r ∂ r + e φ ∂ φ /r + e z ∂ z , velocity u = (u, v), u and v denoting the radial and axial velocity components, and stress tensor σ = −p 1 + η ∇u + (∇u) T , with pressure p and viscosity η. ρ, c p and k are the material density, heat capacity and heat conductivity, respectively, which are all assumed to be constant.
The equations are then made dimensionless applying the following transformation rule, where the subscript '0' denotes the position at z = 0, T max the maximum temperature of the surrounding air (see below) and η min the corresponding viscosity. Note that even though the material does not necessarily reach T max and η min , it still approaches these values closely during the drawing, which makes them well-suited for scaling. α = R 0 /L, the ratio of preform radius to drawing length, is also referred to as the fiber parameter.
At the fiber surface r = R(z), we impose the following boundary conditions, which are also often referred to as the kinematic, free stress and Newton's cooling boundary conditions, respectively. The vector normal to the fiber surface is denoted by n, h denotes the heat transfer coefficient between the fiber material and the surrounding air, assumed to be constant, and T air the temperature profile of the surrounding air. While radiation is neglected in (S3c), we adapt h and T air later by fitting the radius profile of the fiber (see below), which means we are using an effective heat transfer coefficient and air temperature profile to cover both convection and radiation using the simple form of (S3c). Note that it is also possible to derive Eq. (S3c) by linearizing the radiative heat transfer [3].
We assume a centered Gaussian profile for the air temperature in the furnace along the drawing axis, which reads in dimensionless form and is determined by the dimensionless amplitude Λ and width ∆. The domain of modeling, i.e., the drawing length L, is set to 40 cm around the maximum of the ambient temperature T max , which is also used for scaling. With this choice of domain, the temperature at the boundaries is small enough so that the resulting high viscosities lead to negligible fiber deformation.
Cross-section averaged model The variables are expanded in orders of α 2 , i.e., with state vector x = (p, u, η, T ). Analogue to the way Scheid et al. [2] included nonisothermal effects in one-dimensional film casting, we assume a parabolic temperature field in r with averaged temperature T , which fulfills boundary condition (S3c), with Biot number Bi = R 0 h/k. This enables us to cross-section average the leading order of the governing equations (S1) by integration to obtain with Stanton number St = 2 Bi/(αP e), where the Peclet number is given by P e = ρc p v 0 R 0 /k. This derivation also yields v (0) = v (0) (z), i.e., the axial velocity at leading order does not vary along the radius.
Equations (S7) are supplemented by Arrhenius' law, which links the heat equation to the rest of the system and is given by The dimensionless activation energy µ = E a /(R T max ), with dimensional activation energy E a and Gas constantR, quantifies the sensitivity of the viscosity on temperature changes.
The boundary conditions, close the cross-section averaged model.

Higher order terms
We now investigate the α 2 -order of the momentum equation (S1b), which yields To solve this equation for v (1) by numerical integration, we require the condition v (1) = 0, (S11) so as to preserve the fact that at leading order, v (0) is equal to the averaged velocity v .
The α 2 -order of the continuity equation (S1a) then leads to instead of the shear rate, which then yields ∂ t p e = −a η II D m (p e − p c ).
(S14) Figure S1 shows the distribution of √ II D in a PSu fiber. The dimensionless destruction rate a is used as a free fit parameter and is proportional to (v 0 /L) m−1 η min , while m is set to 1.5, following Saphiannikova et al. [4]. The critical concentration p c , below which no conductive pathway exists can be eliminated later.
Integrating Eq. (S14) yields with the radial position being a function of z. Equation (S15) can then be rewritten as The conductivity is then determined using a common percolation law, with initial conductivity of the preform σ in and critical exponent c. Substituting Eq. (S16) in (S17) then leads to for the final conductivity of the fully drawn fiber. The initial conductivity of the undeformed material σ in and the exponent a c are determined by fitting the experimental data, similar to the work of Saphiannikova et al. [4].

Materials
Polysulfone and polycarbonate plates are purchased from Boedecker Plastics, Inc. Pellets of polycarbonate and a polycarbonate/carbon nanotube masterbatch (15 wt%) are purchased from Nanocyl SA and melt-mixed together in a DSM twin-screw micro-compounder at 240 • C to achieve a concentration of 3 wt%; the screw rotation speed is set to 100 rpm and the mixing time is 5 min; the extruded material is then hot pressed into a sheet at 200 • C for 2 minutes with a Lauffer press. The polycarbonate/carbon black nanocomposite is purchased from Westlake Plastics in sheet form (0.25 mm thick). Metallic ribbons (2 mm x 1 mm) of low melting point tin-zinc and bismuth-tin alloys are purchased from Indium Corporation.

Preform fabrication
The polysulfone or polycarbonate plates are cut to the dimensions of the metallic consolidation mold (170 mm x 24 mm) employed and milled with grooves dimensioned to fit the nanocomposite sheets and metallic electrodes. After positioning the materials together, they are thermally welded by hot pressing the preform in the consolidation mold in a Lauffer press. A pressure of 1 N/cm 2 is applied for 1 h at a temperature of 220 • C (185 • C) for preforms with a polysulfone (polycarbonate) cladding.

Thermal drawing
The preforms are thermally drawn into fibers in a custom-made draw tower. The temperatures of the three heating zones of the tube furnace are typically set to 150/260/100 • C for polycarbonate and 200/320/120 • C for polysulfone. The feeding speed is 2 mm/s, and the drawing speed is varied depending on the targeted draw ratio.

Conductivity measurements
The fibers have metallic electrodes, which are in contact with the nanocomposites. The cladding above the metallic electrodes is locally polished away to enable a measurement of the electrical resistance perpendicular to the drawing direction of the nanocomposite layer between two electrodes, using a Keithley 2450 Sourcemeter. From the resistance measurement, the conductivity is then calculated using the geometrical dimensions of the region considered. are obtained from the material specification sheet of the supplier and Comsol Multiphysics material library. The effective heat transfer coefficient h and the effective air temperature profile parameters Λ, ∆ and T max , the latter being necessary for calculating µ, are determined by fitting the predicted fiber radius profile to the experimental observation. In addition to that, we make sure that the experimentally measurable tension necessary to draw the fiber is reproduced by the theoretical prediction. As the fabricated fibers have a rectangular cross-section, we use both the dimensionless width and thickness profiles for fitting. Both profiles exhibit a very similar dimensionless form, which comforts the use of an axisymmetric

PARAMETER DETERMINATION
model. An example of such a fit is given by Fig. S2 and the parameters determined for both cladding materials are given in Table S2.

EVOLUTION OF TEMPERATURE, VISCOSITY AND VELOCITY
We provide here additional information on the axial and radial shape of temperature, viscosity and axial velocity. Figure S3 compares the air temperature profile to the crosssection averaged temperature in a PSu fiber. While the maximum temperature T max is almost reached by the material, justifying the choice of it as reference scale, the material temperature exhibits a delayed reaction to the air temperature, which is especially visible in the cooling part (z > 0.5). The delay is primarily determined by the Stanton number, which is St = 32 here, as it compares the heat transfer between the fiber and the surrounding to the internal convection. A higher St thus means T approaching T air closer.   Figure S4 shows the evolution of cross-section averaged viscosity, axial velocity and strain rate along the drawing axis. Similar to the temperature, the viscosity does not reach its initial value at z = 1, but is six orders of magnitude higher than its lowest value around the middle of the drawing distance. As can be seen in the velocity and strain rate profiles, the fiber deformation is already finished at an earlier position and most of the deformation occurs between z = 0.4 and 0.8.
The insets in Fig. S4 show the radial distribution of the non-averaged quantities at several positions on the drawing axis. While viscosity and strain rate exhibit significant changes along the radius over the entire drawing length, the velocity varies only in a narrow range from the averaged value.