2. The free-boundary problem

The following general set-up closely follows that described by Stokes et al. (Reference Stokes, Buchak, Crowdy and Ebendorff-Heidepriem2014). A preform with cross-sectional area $S_{0}$ containing $M\geqslant 1$ channels is fed vertically down into a hot zone at speed $U_{0}$ and pulled down at a draw speed $U_{1}$ to a cross-sectional area $S_{1}$ over a ‘neck-down length’ $L$ . The shapes of the channels deform during the draw due to the combined effects of surface tension and draw tension. A schematic of the drawing of such a fibre is shown in figure 1. The $x$ -axis is directed vertically down along the axis of the fibre.

As shown in Stokes et al. (Reference Stokes, Buchak, Crowdy and Ebendorff-Heidepriem2014), neglect of inertia is justified because the Reynolds number typical for fibre drawing is very small, $O(10^{-8})$ . We also neglect gravity. Thus, the material making up the preform is assumed to be a viscous Newtonian fluid in the Stokes regime, and its incompressible velocity field $\boldsymbol{u}$ satisfies

(2.1a,b ) $$\begin{eqnarray}-\boldsymbol{{\rm\nabla}}p+{\it\mu}{\rm\nabla}^{2}\boldsymbol{u}=0,\quad \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0,\end{eqnarray}$$

where $p$ is the pressure in the fluid and ${\it\mu}$ is its viscosity, here assumed to be constant and uniform throughout the fluid. The fluid has $M+1$ free surfaces, one for each of the $M$ channels, together with an outer boundary. Using $k=0$ for the outer boundary and $k=1,\dots ,M$ for the $k$ th channel, the dynamic boundary condition on the $k$ th free surface is

(2.2) $$\begin{eqnarray}{\it\sigma}_{ij}n_{j}=-{\it\gamma}{\it\kappa}n_{i}-p_{B}^{(k)}n_{i},\end{eqnarray}$$

where ${\it\kappa}$ is the surface curvature, ${\it\gamma}$ is the surface tension, ${\it\sigma}_{ij}$ is the fluid stress tensor and $p_{B}^{(k)}(t)$ is the pressure on the $k$ th boundary. Without loss of generality, we set $p_{B}^{(0)}=0$ . Under the assumption of quasi-steady flow, the normal velocity $V_{k}$ of each boundary point is given by

(2.3) $$\begin{eqnarray}V_{k}=\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n}.\end{eqnarray}$$

In addition, the boundary conditions at $x=0$ and $x=L$ are $\boldsymbol{u}=U_{0}\boldsymbol{i}$ and $\boldsymbol{u}=U_{1}\boldsymbol{i}$ respectively.

Following Cummings & Howell (Reference Cummings and Howell1999), we introduce a small parameter ${\it\epsilon}$ characterising the slenderness of the fibre,

(2.4) $$\begin{eqnarray}{\it\epsilon}=\frac{\sqrt{S_{0}}}{L}\ll 1,\end{eqnarray}$$

which is found to be well satisfied for fibres typically encountered in practice. We also define the dimensionless surface tension ${\it\gamma}^{\ast }$ in terms of the physical parameters of the problem,

(2.5) $$\begin{eqnarray}{\it\gamma}^{\ast }=\frac{{\it\gamma}L}{{\it\mu}\sqrt{S_{0}}U_{0}}.\end{eqnarray}$$

Cummings & Howell (Reference Cummings and Howell1999) showed that for a slender fibre ( ${\it\epsilon}\ll 1$ ), the three-dimensional problem splits into two simpler problems: a one-dimensional problem in the axial (i.e.  $x$ ) direction (‘stretching’) and a two-dimensional problem in the cross-plane (‘sintering’). The cross-plane problem turns out to be simply the classical two-dimensional problem of surface-tension-driven Stokes flow with unit cross-sectional area, surface tension and viscosity, evolving with respect to a reduced time variable ${\it\tau}$ to be introduced shortly.

The stretching problem is governed by a pair of differential equations in the axial coordinate $x$ and time $t$ . Generalisation of the derivation of Cummings & Howell (Reference Cummings and Howell1999) to allow for the presence of pressurised channels, conservation of mass and momentum in each cross-section requires respectively

(2.6) $$\begin{eqnarray}\displaystyle & S_{t}+(uS)_{x}=0, & \displaystyle\end{eqnarray}$$

(2.7) $$\begin{eqnarray}\displaystyle & \displaystyle 3(Su_{x})_{x}+\frac{1}{2}{\it\gamma}^{\ast }\mathop{\sum }_{k=0}^{M}({\it\Gamma}_{k})_{x}=\mathop{\sum }_{k=0}^{M}\pm (p_{B}^{(k)})_{x}\mathscr{A}_{k}. & \displaystyle\end{eqnarray}$$

$S(x,t)$

$u(x,t)$

$\mathscr{A}_{k}(x,t)$

${\it\Gamma}_{k}(x,t)$

$p_{B}^{(k)}(x,t)$

$k$

$\tilde{x}$

$\tilde{t}$

(2.8) $$\begin{eqnarray}t=\tilde{t},\quad x=\int _{0}^{\tilde{t}}u(\tilde{x},T)\,\text{d}T+\tilde{x},\end{eqnarray}$$

the reduced time ${\it\tau}$ is defined by (Cummings & Howell Reference Cummings and Howell1999; Griffiths & Howell Reference Griffiths and Howell2008)

(2.9) $$\begin{eqnarray}{\it\tau}={\it\gamma}^{\ast }\int _{0}^{\tilde{t}}\frac{\text{d}\tilde{t}}{\sqrt{S}}.\end{eqnarray}$$

In steady fibre drawing, this reduced time ${\it\tau}$ characterising the cross-plane evolution is related to the $x$ coordinate by

(2.10) $$\begin{eqnarray}\frac{\partial x}{\partial {\it\tau}}=\frac{u\sqrt{S}}{{\it\gamma}^{\ast }}.\end{eqnarray}$$

In (2.6)–(2.10), distances in the axial direction have been scaled by $L$ , distances in the cross-plane have been scaled by $\sqrt{S_{0}}$ , the axial velocity has been scaled by $U_{0}$ , and time and pressure have been scaled by $L/U_{0}$ and ${\it\mu}U_{0}/L$ respectively.

In the next section we summarise results from a companion paper (Stokes et al. Reference Stokes, Buchak, Crowdy and Ebendorff-Heidepriem2014), where it is demonstrated that the solution to the axial problem for arbitrary cross-plane geometries can be given in terms of explicit formulae. It should be noted that the formulation of Stokes et al. (Reference Stokes, Buchak, Crowdy and Ebendorff-Heidepriem2014) allows for an axially varying viscosity profile, while here we focus on the case where the fluid viscosity is uniform. Should it be desirable to consider a non-constant viscosity, the required modifications to the model are as in Stokes et al. (Reference Stokes, Buchak, Crowdy and Ebendorff-Heidepriem2014).

3. Explicit formulae for stretching and model predictions

In this section, we focus on the special case in which none of the channels are pressurised relative to the ambient pressure, therefore

(3.1) $$\begin{eqnarray}p_{B}^{(k)}(x,t)=0,\quad k=0,1,\dots ,M.\end{eqnarray}$$

In steady fibre drawing, it is expedient (Stokes et al. Reference Stokes, Buchak, Crowdy and Ebendorff-Heidepriem2014) to introduce an important function

(3.2) $$\begin{eqnarray}H({\it\tau})=\exp \left(-\frac{1}{12}\int _{0}^{{\it\tau}}\tilde{{\it\Gamma}}({\it\tau}^{\prime })\,\text{d}{\it\tau}^{\prime }\right),\end{eqnarray}$$

where $\tilde{{\it\Gamma}}({\it\tau})$ is the total boundary perimeter in the cross-plane, scaled by $\sqrt{S_{0}S({\it\tau})}$ . The total boundary perimeter $\tilde{{\it\Gamma}}({\it\tau})$ can be readily computed as the sum of the perimeters of the individual boundaries once the solution for the cross-plane problem, with total area normalised to be unity, has been found. With $H({\it\tau})$ known from (3.2), for a fluid with constant uniform viscosity the cross-section $S({\it\tau})$ and axial position $x({\it\tau})$ are given explicitly as

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle S({\it\tau})=\left(1-\frac{{\it\sigma}^{\ast }}{{\it\gamma}^{\ast }}\int _{0}^{{\it\tau}}H({\it\tau}^{\prime })\,\text{d}{\it\tau}^{\prime }\right)^{2}\!\!\left.\vphantom{\int }\right/H({\it\tau})^{2}, & \displaystyle\end{eqnarray}$$

(3.4) $$\begin{eqnarray}\displaystyle & \displaystyle x({\it\tau})=-\frac{1}{{\it\sigma}^{\ast }}\log (H({\it\tau})\sqrt{S({\it\tau})}). & \displaystyle\end{eqnarray}$$

${\it\sigma}^{\ast }$

${\it\sigma}$

(3.5) $$\begin{eqnarray}{\it\sigma}=\frac{6{\it\mu}U_{0}S_{0}}{L}{\it\sigma}^{\ast }.\end{eqnarray}$$

Formulae (3.3) and (3.4) give the complete solution, parametrised by ${\it\tau}$ , to the three-dimensional fibre drawing problem from a solution of the cross-plane problem.

It is interesting to point out that, in devising a scheme to apply the predictions of the model described in Fitt et al. (Reference Fitt, Furusawa, Monro, Please and Richardson2002) for the pulling of a capillary tube – that is, a fibre with just one central channel – to the drawing of MOFs with several channels, Kostecki et al. (Reference Kostecki, Ebendorff-Heidepriem, Warren-Smith and Monro2014) were compelled to study the multichannel geometry to find some quantity that best approximates the single internal channel radius appearing in the equations of the Fitt model. They found that, among several other apparently reasonable options, the best choice is to use the total perimeter length of the collection of channels in the MOF cross-plane. Intuitively, this empirical observation is consonant with the fact that the total cross-plane perimeter plays a crucial role in our own formulation just described.

According to the approach advocated by Stokes et al. (Reference Stokes, Buchak, Crowdy and Ebendorff-Heidepriem2014), ${\it\tau}=0$ labels the initial channel configuration, and one chooses the value ${\it\tau}={\it\tau}_{L}$ corresponding to the desired fibre configuration. In this way ${\it\tau}_{L}$ is the duration of the ‘sintering’ process. By ${\it\tau}={\it\tau}_{L}$ , the cross-section must reach the scaled axial distance to the end of the draw, $x=1$ , and the axial fluid velocity must reach the specified scaled draw speed, $u=U_{1}/U_{0}=D$ , where $D$ is the draw ratio. These boundary conditions give two constraints on the four parameters ${\it\tau}_{L}$ , ${\it\gamma}^{\ast }$ , $D$ and ${\it\sigma}^{\ast }$ ,

(3.6a,b ) $$\begin{eqnarray}{\it\sigma}^{\ast }=\log \left(\frac{\sqrt{D}}{H({\it\tau}_{L})}\right),\quad \frac{H({\it\tau}_{L})}{\sqrt{D}}=1-\frac{{\it\sigma}^{\ast }}{{\it\gamma}^{\ast }}\int _{0}^{{\it\tau}_{L}}H({\it\tau}^{\prime })\,\text{d}{\it\tau}^{\prime }.\end{eqnarray}$$

From (2.5) and (3.5), the ratio of the dimensionless parameters ${\it\gamma}^{\ast }$ and ${\it\sigma}^{\ast }$ is related to the ratio of the physical tensions ${\it\gamma}$ and ${\it\sigma}$ by

(3.7) $$\begin{eqnarray}\frac{{\it\sigma}}{{\it\gamma}}=6\sqrt{S_{0}}\frac{{\it\sigma}^{\ast }}{{\it\gamma}^{\ast }}.\end{eqnarray}$$

As discussed in Stokes et al. (Reference Stokes, Buchak, Crowdy and Ebendorff-Heidepriem2014), this is an important result. For drawing of a particular fibre from a given preform, ${\it\tau}_{L}$ and $D$ are known, and since the physical surface tension ${\it\gamma}$ can be considered to be known a priori, (3.6b ) and (3.7) can be solved for the required physical fibre tension ${\it\sigma}$ . Conversely, when $D$ and ${\it\sigma}$ are measured, (3.6b ) and (3.7) can be solved for ${\it\tau}_{L}$ , which gives the final fibre geometry. Thus, use of the measurable and controllable physical fibre tension ${\it\sigma}$ circumvents the need for specific knowledge of the unknown fluid viscosity ${\it\mu}$ and neck-down length $L$ for the draw, both of which can vary with the draw parameters. This was also noted by Chen & Birks (Reference Chen and Birks2013).

If desired, the ratio ${\it\mu}/L$ can be computed by combining (3.5) and (3.6a ). Further, an estimate of the fluid viscosity ${\it\mu}$ can be obtained if, as done by previous authors (Fitt et al. Reference Fitt, Furusawa, Monro, Please and Richardson2002), the neck-down length $L$ is assumed to be equal to the length of the heated zone. This assumption then permits the physical geometry to be determined over the neck-down length. However, we emphasise that it is usually the final geometry that is of interest rather than the geometry throughout the neck-down zone, so that ${\it\gamma}^{\ast }$ and ${\it\sigma}^{\ast }$ need not be computed separately, and there is no need to know or determine the viscosity ${\it\mu}$ or to make any assumptions about the neck-down length $L$ . The slightly more complex case of axially varying viscosity is considered in Stokes et al. (Reference Stokes, Buchak, Crowdy and Ebendorff-Heidepriem2014), leading to the same result for any viscosity profile.

The conditions (3.6a,b ) may be written more succinctly (Stokes et al. Reference Stokes, Buchak, Crowdy and Ebendorff-Heidepriem2014) as the single constraint

(3.8) $$\begin{eqnarray}\frac{1}{\mathscr{Q}}+\frac{\log \mathscr{Q}}{\mathscr{P}}=1,\end{eqnarray}$$

where

(3.9a,b ) $$\begin{eqnarray}\mathscr{P}=\frac{{\it\gamma}^{\ast }}{\displaystyle \int _{0}^{{\it\tau}_{L}}H({\it\tau}^{\prime })\,\text{d}{\it\tau}^{\prime }},\quad \mathscr{Q}=\frac{\sqrt{D}}{H({\it\tau}_{L})}.\end{eqnarray}$$

The choice of ${\it\tau}_{L}$ picks out the final geometry and, in turn, gives the denominators in the definitions (3.9) of $\mathscr{P}$ and $\mathscr{Q}$ . The variable $\mathscr{P}$ gives a measure of the surface tension while $\mathscr{Q}=\exp ({\it\sigma}^{\ast })$ is a measure of the draw tension.

Formula (3.8) represents the balance between draw tension and surface tension required to achieve a desired geometry in a constant-viscosity fluid with no channel pressurisation. A graph of $\mathscr{Q}$ against $\mathscr{P}$ is shown in figure 3. From it we learn that small changes in $\mathscr{P}$ (surface tension) can lead to large corresponding changes in $\mathscr{Q}$ , and hence draw ratio $D$ , for the same desired final configuration. The logarithmic dependence of draw tension ${\it\sigma}^{\ast }$ on $\mathscr{Q}$ means, however, that the corresponding change in draw tension is less sensitive to changes in $\mathscr{P}$ .

Figure 3. The balance of tensions as given by the curve (3.8). The region $\mathscr{Q}>1$ corresponds to positive draw tension. Increasing $\mathscr{P}$ , and hence surface tension so that the cross-plane deformations happen over a shorter time scale, requires a higher $\mathscr{Q}$ (essentially, draw ratio) to achieve the same target geometry.

4. The cross-plane problem

To model the drawing of MOFs of general shape it is necessary to solve for the evolution of a given cross-plane geometry evolving under the effect of surface tension. We now briefly describe the free-boundary cross-plane problem and the formulation of a numerical method that we have devised for its solution. This numerical method is interesting in itself and will be described in detail elsewhere (Buchak & Crowdy Reference Buchak and Crowdy2014). Its main use for the present purposes is simply to test the accuracy of an approximate (elliptical pore) model to be presented in § 5. A boundary integral approach to the same numerical problem has been proposed by Chakravarthy & Chiu (Reference Chakravarthy and Chiu2009), but it differs in several respects from the approach adopted here.

The cross-plane problem for steady fibre drawing is governed by the Stokes equations,

(4.1a,b ) $$\begin{eqnarray}-\boldsymbol{{\rm\nabla}}p+{\it\mu}{\rm\nabla}^{2}\boldsymbol{u}=0,\quad \boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{u}=0,\end{eqnarray}$$

with boundary conditions

(4.2a,b ) $$\begin{eqnarray}{\it\sigma}_{ij}n_{j}=-{\it\gamma}{\it\kappa}n_{i}-p_{B}^{(n)}n_{i},\quad V_{n}=\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{n},\end{eqnarray}$$

where now the flow is two-dimensional and the boundaries are closed curves in the cross-plane. Lengths are scaled by $\sqrt{S_{0}S({\it\tau})}$ , velocities by $\sqrt{S_{0}}U_{0}{\it\gamma}^{\ast }/L$ and pressures by ${\it\mu}U_{0}{\it\gamma}^{\ast }/L\sqrt{S({\it\tau})}$ . The fluid region evolves in the reduced time ${\it\tau}$ defined in (2.9) and has unit area, unit surface tension and unit viscosity (our numerical method can be applied to more general two-dimensional Stokes flow, so for this reason we retain the constants ${\it\mu}$ and ${\it\gamma}$ in the following description).

Two-dimensional incompressible flow can be formulated in terms of a streamfunction ${\it\psi}(x,y)$ , satisfying

(4.3) $$\begin{eqnarray}{\rm\nabla}^{4}{\it\psi}=0,\end{eqnarray}$$

where ${\rm\nabla}^{2}$ denotes the two-dimensional Laplacian operator in the cross-plane. Introducing the complex variable $z$ as the cross-plane coordinate in the usual way, the general solution to (4.3) can be written

(4.4) $$\begin{eqnarray}{\it\psi}=\text{Im}[\overline{z}f(z,{\it\tau})+g(z,{\it\tau})],\end{eqnarray}$$

where $f(z,{\it\tau})$ and $g(z,{\it\tau})$ are two functions analytic in the fluid. The fluid velocity $(u,v)$ , vorticity ${\it\omega}$ , pressure $p$ and rate-of-strain tensor $\unicode[STIX]{x1D626}_{ij}$ are given by (Langlois Reference Langlois1964; Crowdy Reference Crowdy2003a )

(4.5) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{p}{{\it\mu}}-\text{i}{\it\omega}=4f^{\prime }(z,{\it\tau}), & \displaystyle\end{eqnarray}$$

(4.6) $$\begin{eqnarray}\displaystyle & u+\text{i}v=-f(z,{\it\tau})+z\overline{f^{\prime }(z,{\it\tau})}+\overline{g^{\prime }(z,{\it\tau})}, & \displaystyle\end{eqnarray}$$

(4.7) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D626}_{11}+\text{i}\unicode[STIX]{x1D626}_{12}=z\overline{f^{\prime \prime }(z,{\it\tau})}+\overline{g^{\prime \prime }(z,{\it\tau})}. & \displaystyle\end{eqnarray}$$

$n~(n=0,\dots ,M)$

(4.8) $$\begin{eqnarray}\displaystyle & f(z,{\it\tau})+z\overline{f^{\prime }(z,{\it\tau})}+\overline{g^{\prime }(z,{\it\tau})}=-{\textstyle \frac{1}{2}}\text{i}{\it\gamma}z_{s}+{\textstyle \frac{1}{2}}p_{B}^{(n)}z+A_{n}({\it\tau}), & \displaystyle\end{eqnarray}$$

(4.9) $$\begin{eqnarray}\displaystyle & \displaystyle \text{Re}\left\{\frac{\text{d}z}{\text{d}{\it\tau}}\text{i}\overline{z}_{s}\right\}=\text{Re}\left\{\left(-f(z,{\it\tau})+z\overline{f^{\prime }(z,{\it\tau})}+\overline{g^{\prime }(z,{\it\tau})}\right)\text{i}\overline{z}_{s}\right\}, & \displaystyle\end{eqnarray}$$

$z_{s}=\text{d}z/\text{d}s$

$s$

$z(s)$

$\{A_{n}({\it\tau})\,|\,n=0,\dots ,M\}$

The system (4.1) and (4.2) is a highly nonlinear two-dimensional free-boundary problem. Remarkably, several exact solutions to it are known for special geometries (Hopper Reference Hopper1990; Richardson Reference Richardson1992; Tanveer & Vasconcelos Reference Tanveer and Vasconcelos1995; Cummings, Howison & King Reference Cummings, Howison and King1999; Crowdy Reference Crowdy2003a ,Reference Crowdy b ). Numerical methods must be used, however, for general initial conditions. One option is to make use of standard boundary integral methods (Pozrikidis Reference Pozrikidis1992; Chakravarthy & Chiu Reference Chakravarthy and Chiu2009), but we eschew this approach in favour of one based on a conformal mapping representation of the boundaries. One reason is that the boundary integral formulation for compressible bubbles, with changing area, suffers from certain technical difficulties involving the appearance of hypersingular integrals in the standard formulation (Pozrikidis Reference Pozrikidis2001, Reference Pozrikidis2003); a second reason is that it is precisely such conformal mapping formulations that give rise to the known exact solutions mentioned above, in particular those for compressible bubbles given by Crowdy (Reference Crowdy2003a ) and which lie at the heart of the model to be described in § 5.

4.1. Numerical method

To track the fluid boundaries in the cross-plane we introduce a ${\it\tau}$ -dependent conformal mapping $z({\it\zeta},{\it\tau})$ from a canonical circular region in a complex ${\it\zeta}$ -plane to the evolving fluid region. This canonical preimage region is inside the circle $|{\it\zeta}-{\it\delta}_{0}({\it\tau})|=q_{0}({\it\tau})$ but outside the circles $|{\it\zeta}-{\it\delta}_{n}({\it\tau})|=q_{n}({\it\tau})~(n=1,\dots ,M)$ , where the circle centres ${\it\delta}_{n}({\it\tau})\in \mathbb{C}$ and radii $q_{n}({\it\tau})\in \mathbb{R}$ evolve in time. For convenience, we take ${\it\delta}_{0}({\it\tau})=0$ and $q_{0}({\it\tau})=1$ . The conformal mapping can be represented in general as

(4.10) $$\begin{eqnarray}z({\it\zeta},{\it\tau})=\mathop{\sum }_{n=0}^{M}\mathop{\sum }_{k}a_{k}^{(n)}({\it\tau})\left(\frac{{\it\zeta}-{\it\delta}_{n}({\it\tau})}{q_{n}({\it\tau})}\right)^{k},\end{eqnarray}$$

where for the outer boundary ( $n=0$ ) we take $k=0,\dots ,\infty$ , while for the inner boundaries ( $n=1,\dots ,M$ ) we take $k=-1,\dots ,-\infty$ . The coefficients $\{a_{k}^{(n)}({\it\tau})\}$ as well as the parameters $\{{\it\delta}_{n}({\it\tau}),q_{n}({\it\tau})\,|\,n=1,\dots ,M\}$ characterising the preimage circles must be found. With the introduction of this conformal mapping we can now define the composed functions

(4.11a,b ) $$\begin{eqnarray}F({\it\zeta},{\it\tau})\equiv f(z({\it\zeta},{\it\tau}),{\it\tau}),\quad G({\it\zeta},{\it\tau})\equiv g^{\prime }(z({\it\zeta},{\it\tau}),{\it\tau}),\end{eqnarray}$$

which are determined, at each ${\it\tau}$ value, by the stress boundary condition (4.8) provided that some additional degrees of freedom are chosen correctly (Buchak & Crowdy Reference Buchak and Crowdy2014). Once $F({\it\zeta},{\it\tau})$ and $G({\it\zeta},{\it\tau})$ are known, the normal velocity of the boundary can be computed from the kinematic condition (4.9) and the boundary position updated in ${\it\tau}$ using a fourth-order Runge–Kutta scheme. A detailed explanation and performance evaluation of this numerical spectral method for computing the cross-plane evolution will be given elsewhere (Buchak & Crowdy Reference Buchak and Crowdy2014).

To validate the numerical method, we checked it against known analytical solutions for the two-dimensional sintering of $N$ equal circular cylinders arranged in a doubly connected annular arrangement (Crowdy Reference Crowdy2003b ). We also verified that our method successfully captures the evolution of a thin viscida of fluid, as predicted by a set of asymptotic equations derived by Griffiths & Howell (Reference Griffiths and Howell2007).

The numerical scheme just described can, in principle, be used to solve the cross-plane problem for any given initial geometry. Here, however, we use it only to test the validity and accuracy of the model introduced in the next section.

5. The elliptical pore model

How can we find initial geometries that, under suitable draw conditions, will produce an MOF with the desired target geometry? Previous authors (Yarin Reference Yarin1995; Griffiths & Howell Reference Griffiths and Howell2007, Reference Griffiths and Howell2008) and the present authors (Stokes et al. Reference Stokes, Buchak, Crowdy and Ebendorff-Heidepriem2014) have addressed this question for special geometries, but it is desirable to be able to handle more general geometries relevant to MOF applications.

Since this is clearly an inverse problem, one option is to consider running the two-dimensional free-surface Stokes flow problem backwards in time, starting with the desired configuration, with the hope of backing out an appropriate initial geometry. However, a little consideration reveals that this inverse problem is ill-posed: surface tension smoothes out ripples in an interface so, in reverse time, any numerical method will be subject to contamination by the growth of initially small numerical inaccuracies. The authors have given a concrete illustrative example of this ill-posedness (Stokes et al. Reference Stokes, Buchak, Crowdy and Ebendorff-Heidepriem2014); Yarin (Reference Yarin1995) has remarked that ‘computers are practically useless in solving such problems’. Griffiths & Howell (Reference Griffiths and Howell2008, Reference Griffiths and Howell2009) have acknowledged the difficulty of a backwards-time calculation. It is therefore crucial, in order to solve this inverse problem in any meaningful way, to strategically constrain the class of possible solutions by adding additional information on the type of initial condition one seeks.

Yarin (Reference Yarin1995) found that truncation of a Fourier series provided the necessary constraint for solving the inverse problems there considered. Griffiths & Howell (Reference Griffiths and Howell2008, Reference Griffiths and Howell2009), in their consideration of the drawing of thin-walled non-axisymmetric tubes, presented an asymptotic model of the evolution of their cross-plane geometries which allows the solution to the cross-plane problem to be written down in analytical form by, in essence, tracking only the centreline of a thin viscida. This approach has the subsidiary – and highly desirable – advantage that time can be reversed in a stable manner and predictions made as to initial conditions giving rise to required end-state geometries. In our own prior work (Stokes et al. Reference Stokes, Buchak, Crowdy and Ebendorff-Heidepriem2014) we made use of various known exact solutions of the cross-plane problem for axisymmetric and non-axisymmetric tubes (of arbitrary wall thickness) to isolate solutions of this inverse problem.

For the geometries typically associated with multichannel MOFs, the models described above cannot be applied. Therefore, we now present a different asymptotic model of the cross-plane flow which we have found to be suited to broad classes of multichannel MOF geometries. The model affords us similar advantages to the models of Griffiths & Howell (Reference Griffiths and Howell2008) and Stokes et al. (Reference Stokes, Buchak, Crowdy and Ebendorff-Heidepriem2014) in that it admits an analytic formulation (thereby obviating the need for a full numerical simulation of the kind expounded in § 4) and it can be run backwards in time in a stable manner to give predictions for initial conditions giving rise to the required MOF geometries at the end of the draw.

The main idea is to employ a so-called elliptical pore model (henceforth, EPM) of two-dimensional interacting bubbles in two-dimensional Stokes flows expounded earlier, and in a different context, by one of the authors (Crowdy Reference Crowdy2004). The basis of the EPM rests on the mathematical observation (Crowdy Reference Crowdy2003a ) that an isolated compressible elliptical bubble (or ‘pore’) in a general time-dependent linear Stokes flow remains precisely elliptical in shape under evolution according to the Stokes equations. Given this surprising fact, it was proposed (Crowdy Reference Crowdy2004) to model multibubble interactions – in the spirit of inner–outer matched asymptotic expansions – by assuming that, in a far-field sense, each compressible bubble which is generally shrinking under the effects of surface tension will be well modelled as an effective point sink centred at its centroid. Hence, provided that all bubbles are sufficiently far apart (relative to a typical length scale of the bubbles), any given bubble will evolve as an elliptical bubble in the time-dependent irrotational linear strain produced by the superposition of point sinks representing the flow induced by its neighbours. Numerical evidence given by Crowdy (Reference Crowdy2004) confirmed that this model gives good agreement with the full dynamical evolution of certain test cases.

5.1. Generalised EPM

The original EPM of Crowdy (Reference Crowdy2004) proposed to treat the pores, at long-range distances, as pure point sinks in an unbounded region. For MOFs we need to extend this model in two ways. First, to more accurately model the channel interactions, we include contributions from another type of Stokes flow singularity – the stresslet (Pozrikidis Reference Pozrikidis1992) – which we describe below. Second, the bounded geometry of an MOF cross-section requires us to incorporate the effects of an outer boundary in our model.

The components of our model are sketched in figure 4. The cross-section of an MOF is characterised by $M$ elliptical channels with centroids located at $\mathscr{Z}_{n}({\it\tau})\in \mathbb{C}~(n=1\dots \,M)$ in the cross-plane. We assume that the outer boundary begins and stays circular and centred at the origin as it evolves – an assumption consistent with many MOF geometries used in practice which often comprise rotationally symmetric arrays of channels. Hence, the outer boundary can be described by

(5.1) $$\begin{eqnarray}|z|=R({\it\tau}),\end{eqnarray}$$

for some radius $R({\it\tau})$ to be determined.

Figure 4. Schematic illustrating the idea of the generalised EPM. Each channel, labelled by $n$ , in the fibre cross-section is modelled as an elliptical bubble in an ambient flow characterised by a local strain rate $k_{n}$ , vorticity ${\it\omega}_{n}$ and pressure $p_{n}$ (not shown). This local flow is induced by the other channels and the outer boundary, assumed to be well separated from channel $n$ .

Each channel is modelled, in its far field, as the combination of a source of strength $m_{n}({\it\tau})$ and a stresslet of strength ${\it\lambda}_{n}({\it\tau})$ ; note that it is easily shown that the net force on the fluid due to the presence of the channel is zero, thereby precluding any use of Stokeslet singularities to model them. In terms of the Goursat functions introduced in (4.4), a point source/sink of strength $m_{n}$ at $\mathscr{Z}_{n}$ is represented as

(5.2a,b ) $$\begin{eqnarray}f(z,{\it\tau})=0,\quad g^{\prime }(z,{\it\tau})=\frac{m_{n}}{2{\rm\pi}(z-\mathscr{Z}_{n})},\end{eqnarray}$$

which, from formula (4.6), gives the velocity field

(5.3) $$\begin{eqnarray}u+\text{i}v=\frac{m_{n}}{2{\rm\pi}\overline{(z-\mathscr{Z}_{n})}}.\end{eqnarray}$$

A point stresslet of strength ${\it\lambda}_{n}$ is represented as

(5.4a,b ) $$\begin{eqnarray}f(z,{\it\tau})=\frac{{\it\lambda}_{n}}{2{\rm\pi}(z-\mathscr{Z}_{n})},\quad g^{\prime }(z,{\it\tau})=\frac{{\it\lambda}_{n}\overline{\mathscr{Z}_{n}}}{2{\rm\pi}(z-\mathscr{Z}_{n})^{2}},\end{eqnarray}$$

which produces the velocity field

(5.5) $$\begin{eqnarray}u+\text{i}v=-\frac{{\it\lambda}_{n}}{2{\rm\pi}(z-\mathscr{Z}_{n})}-\frac{\overline{{\it\lambda}_{n}}(z-\mathscr{Z}_{n})}{2{\rm\pi}\overline{(z-\mathscr{Z}_{n})}^{2}}.\end{eqnarray}$$

We have not found this complex variable form of the stresslet singularity documented in the standard textbook literature, but it has been used previously by one of the authors (Crowdy & Or Reference Crowdy and Or2010) in modelling low-Reynolds-number swimmers near walls. Together, $f(z,{\it\tau})$ and $g^{\prime }(z,{\it\tau})$ give a velocity field that decays like $|z|^{-1}$ in the far field. Since the channels generally shrink in size, we expect $m_{n}<0$ . The stresslet strengths ${\it\lambda}_{n}({\it\tau})$ are generally complex quantities (reflecting an amplitude and orientation).

The scalings for lengths, velocities and pressures are the same as in § 4. To keep the formulation of the model as general as possible, we allow each channel to be pressurised with a given pressure $p_{B}^{(n)}({\it\tau})$ , and use $p_{B}^{(0)}({\it\tau})$ to denote the ambient pressure exterior to the fibre. For the specific purposes of the present paper, however, we later choose $p_{B}^{(n)}({\it\tau})=0$ for $n=0,\dots ,M$ (pressurisation effects are studied separately in a forthcoming paper (Chen et al. Reference Chen, Stokes, Buchak, Crowdy and Ebendorff-Heidepriem2015)).

5.2. Global flow

The EPM assumes the large-scale description of the flow to be given by (4.4) with

(5.6) $$\begin{eqnarray}\displaystyle & \displaystyle f(z,{\it\tau})=\underbrace{\mathop{\sum }_{n=1}^{M}\frac{{\it\lambda}_{n}({\it\tau})}{2{\rm\pi}(z-\mathscr{Z}_{n}({\it\tau}))}}_{\text{stresslets}}+\underbrace{\frac{1}{4}P({\it\tau})z}_{\text{outer}\;\text{boundary}\,\text{effect}}, & \displaystyle\end{eqnarray}$$

(5.7) $$\begin{eqnarray}\displaystyle & \displaystyle g^{\prime }(z,{\it\tau})=\underbrace{\mathop{\sum }_{n=1}^{M}\frac{m_{n}({\it\tau})}{2{\rm\pi}(z-\mathscr{Z}_{n}({\it\tau}))}}_{\text{sinks}}+\underbrace{\mathop{\sum }_{n=1}^{M}\frac{{\it\lambda}_{n}({\it\tau})\overline{\mathscr{Z}_{n}({\it\tau})}}{2{\rm\pi}(z-\mathscr{Z}_{n}({\it\tau}))^{2}}}_{\text{stresslets}}. & \displaystyle\end{eqnarray}$$

The values of $m_{n}$ and ${\it\lambda}_{n}$ are determined below by considering the evolution of each individual channel. The original EPM formulation of Crowdy (Reference Crowdy2004) contained only the sink terms labelled in (5.7), but inclusion of the additional stresslet terms in (5.6) and (5.7) means that all physically admissible singularities with velocities decaying like $|z|^{-1}$ are now incorporated into the model. Because the cross-section of an MOF is a bounded region, the term for the uniform pressure in (5.6) is required to account for the influence of the outer boundary on the inner channels.

5.3. Outer boundary evolution

For the flow described by (5.6) and (5.7), the dynamic (4.8) and kinematic (4.9) boundary conditions on the outer boundary give respectively

(5.8) $$\begin{eqnarray}\displaystyle & \displaystyle P({\it\tau})-p_{B}^{(0)}+\frac{\mathscr{M}}{{\rm\pi}R^{2}}=\frac{1}{R}, & \displaystyle\end{eqnarray}$$

(5.9) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}R}{\text{d}{\it\tau}}=\frac{\mathscr{M}}{2{\rm\pi}R}, & \displaystyle\end{eqnarray}$$

where

(5.10) $$\begin{eqnarray}\mathscr{M}=\mathop{\sum }_{n=1}^{M}m_{n}\end{eqnarray}$$

is the combined sink strength of the channels and where, for consistency with the assumption of a purely circular outer boundary, we need

(5.11) $$\begin{eqnarray}\frac{1}{R}\mathop{\sum }_{n=1}^{M}{\it\lambda}_{n}=0.\end{eqnarray}$$

The latter condition is certainly expected to be satisfied by MOFs of practical relevance, involving rotationally symmetric arrays of channels (the modelling of more complicated outer boundary shapes is reserved as a topic for further investigation). Equation (5.8) relates a ‘bulk pressure’ $P({\it\tau})$ to the ambient pressure $p_{B}^{(0)}({\it\tau})$ , while (5.9) gives the evolution of the outer boundary. These equations are equivalent to modelling the effect of the channels on the outer boundary as a point source/sink of strength $\mathscr{M}({\it\tau})$ situated at the centre of the outer boundary circle. In turn, the internal channels feel the effect of the outer boundary through a time varying bulk pressure $P({\it\tau})$ coupled to the outer boundary motion through (5.8).

Equations (5.8)–(5.11) are derived in detail in appendix A.

5.4. Evolution of channels

To determine how the elliptical channels evolve, we must resolve local details of the flow. Following Crowdy (Reference Crowdy2004), the evolution of channel $n$ can be approximated by studying the local field near $\mathscr{Z}_{n}$ . Taylor expansion of (5.6) and (5.7) produces

(5.12) $$\begin{eqnarray}\left.\begin{array}{@{}rcl@{}}f(z) & = & \displaystyle \left(\frac{1}{4}P\mathscr{Z}_{n}+\mathop{\sum }_{j\neq n}\frac{{\it\lambda}_{j}}{2{\rm\pi}(\mathscr{Z}_{n}-\mathscr{Z}_{j})}\right)+\left(\frac{1}{4}P-\mathop{\sum }_{j\neq n}\frac{{\it\lambda}_{j}}{2{\rm\pi}(\mathscr{Z}_{n}-\mathscr{Z}_{j})^{2}}\right)\left(z-\mathscr{Z}_{n}\right)\\ & & +\,\displaystyle \left(\mathop{\sum }_{j\neq n}\frac{{\it\lambda}_{j}}{2{\rm\pi}(\mathscr{Z}_{n}-\mathscr{Z}_{j})^{3}}\right)\left(z-\mathscr{Z}_{n}\right)^{2}+\cdots \,,\\ g^{\prime }(z) & = & \displaystyle \left(\mathop{\sum }_{j\neq n}\frac{m_{j}}{2{\rm\pi}(\mathscr{Z}_{n}-\mathscr{Z}_{j})}+\mathop{\sum }_{j\neq n}\frac{{\it\lambda}_{j}\overline{\mathscr{Z}_{j}}}{2{\rm\pi}(\mathscr{Z}_{n}-\mathscr{Z}_{j})^{2}}\right)\\ & & +\,\displaystyle \left(-\mathop{\sum }_{j\neq n}\frac{m_{j}}{2{\rm\pi}(\mathscr{Z}_{n}-\mathscr{Z}_{j})^{2}}-\mathop{\sum }_{j\neq n}\frac{{\it\lambda}_{j}\overline{\mathscr{Z}_{j}}}{{\rm\pi}(\mathscr{Z}_{n}-\mathscr{Z}_{j})^{3}}\right)(z-\mathscr{Z}_{n})+\cdots \,.\end{array}\right\}\end{eqnarray}$$

A local expansion of the complex velocity $u+\text{i}v$ about the centre $\mathscr{Z}_{n}$ of the $n$ th channel can then be determined from (4.6). The velocity of the channel centroid is given by the constant term in this expansion,

(5.13) $$\begin{eqnarray}\overline{\frac{\text{d}\mathscr{Z}_{n}}{\text{d}{\it\tau}}}=\overline{u(\mathscr{Z}_{n})}=-\mathop{\sum }_{j\neq n}\frac{{\it\lambda}_{j}\overline{(\mathscr{Z}_{n}-\mathscr{Z}_{j})}}{2{\rm\pi}(\mathscr{Z}_{n}-\mathscr{Z}_{j})^{2}}+\mathop{\sum }_{j\neq n}\frac{m_{j}}{2{\rm\pi}(\mathscr{Z}_{n}-\mathscr{Z}_{j})}-\overline{\mathop{\sum }_{j\neq n}\frac{{\it\lambda}_{j}}{2{\rm\pi}(\mathscr{Z}_{n}-\mathscr{Z}_{j})}}.\end{eqnarray}$$

The linear terms in the expansions of $f(z)$ and $g^{\prime }(z)$ , with the origin moved to $\mathscr{Z}_{n}$ , determine how the $n$ th channel shape evolves. A simple calculation leads to

(5.14) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{1}{4}(p_{n}-\text{i}{\it\omega}_{n})=\frac{1}{4}P-\mathop{\sum }_{j\neq n}\frac{{\it\lambda}_{j}}{2{\rm\pi}(\mathscr{Z}_{n}-\mathscr{Z}_{j})^{2}}, & \displaystyle\end{eqnarray}$$

(5.15) $$\begin{eqnarray}\displaystyle & \displaystyle k_{n}=-\mathop{\sum }_{j\neq n}\frac{m_{j}}{2{\rm\pi}(\mathscr{Z}_{n}-\mathscr{Z}_{j})^{2}}-\mathop{\sum }_{j\neq n}\frac{{\it\lambda}_{j}\overline{(\mathscr{Z}_{j}-\mathscr{Z}_{n})}}{{\rm\pi}(\mathscr{Z}_{n}-\mathscr{Z}_{j})^{3}}, & \displaystyle\end{eqnarray}$$

giving the ambient pressure, vorticity and linear strain rate ‘seen’ by the $n$ th channel in its far field. It should be noted that all quantities (5.13)–(5.15) are invariant under a change of origin.

To close the model, we must determine how the shape of the $n$ th channel evolves in response to the far-field conditions (5.14) and (5.15). Following Crowdy (Reference Crowdy2004), the elliptical cross-section of the $n$ th channel is described by a time-dependent conformal mapping $z_{n}({\it\zeta},{\it\tau})$ of a parametric circle $|{\it\zeta}|=1$ . Each has the form

(5.16) $$\begin{eqnarray}z_{n}({\it\zeta},{\it\tau})=\mathscr{Z}_{n}({\it\tau})+\frac{{\it\alpha}_{n}({\it\tau})}{{\it\zeta}}+{\it\beta}_{n}({\it\tau}){\it\zeta},\end{eqnarray}$$

where $\mathscr{Z}_{n}({\it\tau})$ is the centroid position while ${\it\alpha}_{n}({\it\tau})\in \mathbb{R}$ and ${\it\beta}_{n}({\it\tau})\in \mathbb{C}$ are parameters encoding its orientation, area and eccentricity, with $|{\it\beta}_{n}|<{\it\alpha}_{n}$ . It is easy to check that (5.16) transplants the interior of the unit ${\it\zeta}$ -disc to the exterior of an ellipse. Indeed, if $a_{n}$ and $b_{n}$ are the semi-major and semi-minor axes of the elliptical cross-section of the channel, then

(5.17a,b ) $$\begin{eqnarray}a_{n}={\it\alpha}_{n}+|{\it\beta}_{n}|,\quad b_{n}={\it\alpha}_{n}-|{\it\beta}_{n}|,\end{eqnarray}$$

while its area is

(5.18) $$\begin{eqnarray}{\rm\pi}a_{n}b_{n}={\rm\pi}({\it\alpha}_{n}^{2}-|{\it\beta}_{n}|^{2}).\end{eqnarray}$$

The orientation angle of the semi-major axis of the ellipse to the horizontal axis in the cross-plane is $(1/2)\arg [{\it\beta}_{n}]$ .

The evolution of a compressible elliptical bubble in an ambient linear flow was found in analytical form by Crowdy (Reference Crowdy2003a ) (see also Pozrikidis (Reference Pozrikidis2003) for a purely numerical treatment of the same problem). Those results can be applied to the cross-section of an elliptical channel with far-field conditions on $f(z)$ and $g(z)$ taken to be

(5.19) $$\begin{eqnarray}\displaystyle & \displaystyle f(z)\sim \frac{1}{4}(p_{n}-p_{B}^{(n)}-\text{i}{\it\omega}_{n})z+\frac{{\it\lambda}_{n}}{2{\rm\pi}z}+\cdots \,, & \displaystyle\end{eqnarray}$$

(5.20) $$\begin{eqnarray}\displaystyle & \displaystyle g^{\prime }(z)\sim k_{n}z+\frac{m_{n}}{2{\rm\pi}z}+\cdots \,, & \displaystyle\end{eqnarray}$$

where the coefficients $p_{n}-\text{i}{\it\omega}_{n}$ and $k_{n}$ are read off from (5.14) and (5.15). The dynamic and kinematic boundary conditions (4.8) and (4.9) determine the evolution of the parameters ${\it\alpha}_{n}({\it\tau})$ and ${\it\beta}_{n}({\it\tau})$ (Crowdy Reference Crowdy2003a ),

(5.21) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\alpha}_{n}}{\partial {\it\tau}}=-{\it\alpha}_{n}I_{n}(0)-\frac{1}{2}{\it\alpha}_{n}(p_{n}-p_{B}^{(n)}), & \displaystyle\end{eqnarray}$$

(5.22) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\beta}_{n}}{\partial {\it\tau}}=-{\it\beta}_{n}I_{n}(0)+\frac{1}{2}{\it\beta}_{n}(p_{n}-p_{B}^{(n)})+2\overline{k}_{n}{\it\alpha}_{n}+\text{i}{\it\omega}_{n}{\it\beta}_{n}, & \displaystyle\end{eqnarray}$$

with $m_{n}({\it\tau})$ and ${\it\lambda}_{n}({\it\tau})$ for each channel related to ${\it\alpha}_{n}({\it\tau})$ and ${\it\beta}_{n}({\it\tau})$ by

(5.23) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{m_{n}}{{\rm\pi}}=\frac{\text{d}}{\text{d}{\it\tau}}({\it\alpha}_{n}^{2}-|{\it\beta}_{n}|^{2}), & \displaystyle\end{eqnarray}$$

(5.24) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{{\it\lambda}_{n}}{{\rm\pi}}={\it\beta}_{n}\frac{\partial {\it\alpha}_{n}}{\partial {\it\tau}}-{\it\alpha}_{n}\frac{\partial {\it\beta}_{n}}{\partial {\it\tau}}+\text{i}{\it\alpha}_{n}{\it\beta}_{n}{\it\omega}_{n}+2{\it\alpha}_{n}{\it\beta}_{n}I_{n}(0)-\frac{1}{2}{\it\alpha}_{n}^{2}I_{n}^{\prime \prime }(0). & \displaystyle\end{eqnarray}$$

Expressions for $I_{n}(0)$ and $I_{n}^{\prime \prime }(0)$ are given in appendix B. It is apparent from (5.23) and (5.24) that the coefficient ${\it\lambda}_{n}$ may be complex but $m_{n}$ must be real, showing that the channel acts like a superposed stresslet and a source/sink, but has no rotlet contribution.

It is natural to ask what the effect would be of expanding $f(z,t)$ and $g^{\prime }(z,t)$ in (5.12), and hence the velocity field, beyond the linear terms included above. In that case, it is no longer true that an initially elliptical channel sitting in such a higher-order, now nonlinear, ambient flow will remain elliptical. However, some important mathematical generalisations can still be made, and this is discussed further in § 10.

It is instructive to point out that, for a single circular channel at the centre of a circular fibre (i.e. a concentric annular tube), the equations just described retrieve exactly the equations for a pressurised annular fibre that one would obtain by a direct analysis of the equations of motion with radial symmetry (that is, without proceeding via a complex variable and conformal mapping formulation). In other words, the EPM is exact in this case, not an approximation. This is shown in detail in appendix C. The EPM can therefore be understood as a generalised system for the case where the cross-plane comprises more than one channel. It is an approximate system in that case, but one that gives excellent agreement with the exact dynamics provided that the channels are well separated and sufficiently far from the other boundary, as we survey in the next section.

6. Effectiveness of the EPM

The model just described leads to a reduction of the full nonlinear free-boundary problem to a coupled system of ordinary differential equations governing the evolution of the parameters $\mathscr{Z}_{n}({\it\tau})$ , ${\it\alpha}_{n}({\it\tau})$ , ${\it\beta}_{n}({\it\tau})$ and $R({\it\tau})$ . A linear system of equations must be solved at each time step in order to calculate the time derivatives of the parameters. Appendix B illustrates how this can be implemented. We have also made available freely downloadable MATLAB scripts (http://wwwf.imperial.ac.uk/∼dgcrowdy/) for general use; these scripts carry out the prescription given in appendix B.

Henceforth we impose zero channel pressurisation (3.1). As a first check we verified that the model agreed with a known exact solution for a concentric circular geometry (equations (3.2) of van de Vorst (Reference van de Vorst1993), with $R_{i}=0.5$ and $R_{o}=1$ ).

The numerical method of § 4 can be used to test the EPM more generally. We take as a first test configuration the geometry sketched in figure 5, comprising $M$ elliptical channels characterised by the ratio $r=|\mathscr{Z}_{n}|/R$ specifying how close the channels are to the outer boundary. In figure 6, we compare our model with the full solution for a hypothetical MOF whose preform is given by this geometry. (The initial map from a circular domain to the initial configuration is obtained using an extension of the method of Fornberg (Reference Fornberg1980) developed by Kropf (Reference Kropf2009).) Although surface tension causes the shape and size of the channels to change, it is apparent that they remain nearly elliptical in shape. It is also evident that the EPM accurately predicts the time evolution of these ellipses.

Figure 5. The configuration used to test the EPM in figures 6–8. The geometry is characterised by $M$ , the number of channels, and $r=|\mathscr{Z}_{n}|/R$ , the ratio of the radius of the channel centroids to the radius of the outer boundary. Here, $R$ is chosen so that the fluid region has unit area.

Figure 6. Snapshots of the initial configuration (a) and two later configurations (b,c) of an MOF whose initial geometry is shown in figure 5, with $M=4$ and $r=0.4$ : (a)  ${\it\tau}=0$ , (b)  ${\it\tau}=0.040$ , (c)  ${\it\tau}=0.080$ . The solid lines show the evolution according to the EPM; the superposed dashed lines are the full solution using the method of § 4. The middle and right-hand columns show close-ups at each time.

Quantitative measures of the accuracy of the EPM are shown in figures 7, 8, where we plot various quantities associated with one of the channels for the test configuration. The agreement is best when the channels are not too close to each other, or to the outer boundary, which is to be expected given the nature of our approximations.

Figure 7. Testing the accuracy of the model. Instantaneous rates of change of the centroid location (more precisely, its modulus), area $A$ , ratio ${\it\epsilon}$ of major to minor axes of a channel and radius of the outer boundary for the MOF illustrated in figure 5, as the number of channels $M$ is varied, for $r=|\mathscr{Z}_{n}|/R=0.4$ . The solid lines show the value of each quantity according to the EPM; the dashed lines show the value obtained from the full numerical solution. The agreement between the EPM and the full solution is best when $M$ is moderate, that is, when the channels are not too close to each other.

Figure 8. Instantaneous rates of change of the quantities plotted in figure 7 as $r=|\mathscr{Z}_{n}|/R$ is varied, for $M=4$ channels: (a) centroid location, (b) area, (c) axis ratio, (d) outer boundary radius. The solid lines show the value of each quantity according to the EPM; the dashed lines show the value obtained from the full numerical solution. The agreement between the EPM and the full solution is best when $r$ is small, that is, when the channels are far from the outer boundary.

In figure 9, we illustrate the accuracy of the EPM for a hexagonal array, a more realistic geometry for an MOF. In this example, differential straining in the cross-plane leads to certain channels becoming more eccentric than others. In particular, the channels along the innermost ring tend to become extended along the azimuthal direction, while those at the corners of the hexagon on the outermost ring become extended radially. The channels remain roughly elliptical under evolution, and their movement and deformation are captured well by the EPM. This figure shows that the agreement between the model and the full solution is good for a realistic MOF.

Figure 9. Snapshots of the time evolution of an MOF obtained using the EPM (solid lines) and using the numerical method (dashed lines): (a)  ${\it\tau}=0$ , (b)  ${\it\tau}=0.025$ , (c)  ${\it\tau}=0.050$ . The middle and right-hand columns show close-ups at each time.

7. Preliminary experimental validation of the EPM

The authors are currently engaged in a suite of experimental fabrication runs aimed at testing the scope of the mathematical models presented here and in our related modelling work (Stokes et al. Reference Stokes, Buchak, Crowdy and Ebendorff-Heidepriem2014; Chen et al. Reference Chen, Stokes, Buchak, Crowdy and Ebendorff-Heidepriem2015). Full details of those results will be presented elsewhere, but it is appropriate to include here preliminary experimental evidence that the EPM (and future variants thereof) has the potential to faithfully predict the kind of qualitative deformations seen in real fibredraws.

For our tests, we chose a preform consisting of six circular channels arranged in a triangle with a circular outer boundary. This configuration is shown in figure 10(a), along with the predicted configurations at ${\it\tau}=0.025,0.050$ and 0.075 in figure 10(b–d). This same initial configuration was used for an experiment, with six holes of diameter 2.8 mm being drilled into an F2 glass cylinder with an outer diameter of 30 mm. The resulting preform was then drawn into fibre using a fixed feed speed and various tensions and draw speeds. Photographs of the cross-section of the fibre were taken for each set of draw parameters. Two photographs, corresponding to two different choices of the draw parameters, are shown in figure 10(e,f).

Figure 10. Comparison of the EPM with experiment. (a–d) A circular preform with six circular channels (a), and its evolution predicted by the EPM at different times: (a)  ${\it\tau}=0$ , (b)  ${\it\tau}=0.025$ , (c)  ${\it\tau}=0.050$ , (d)  ${\it\tau}=0.075$ . These configurations correspond to possible fibres drawn from this preform. (e,f) Photographs of the fibre cross-section from an experimental draw of a preform with the configuration shown in (a), for two different values of the draw speed and fibre tension. The feed speed was held constant. The scale bars show $50~{\rm\mu}\text{m}$ .

The qualitative agreement is generally good, although it should be noted that, in the second experimental cross-section, the channels close to the centre are larger than those nearer the edge. This may be due to self-pressurisation, perhaps compounded by a radial temperature gradient; neither of these possibilities has been incorporated into our model. Although not a quantitative test, the agreement shows that the overall channel deformation characteristics predicted by the model are also obtained experimentally in fibre drawing, giving confidence in our modelling approach. Quantitative comparisons between the model and experiment are part of ongoing work to be reported in the future.

8. The scope of the EPM

The EPM is only technically valid when the channels are well separated and sufficiently far from the outer boundary, but we have found that its predictions compare favourably with a full numerical solution even for channels that are so close together that the underlying asymptotic separation of ‘inner’ and ‘outer’ scales is not expected to pertain. However, in such cases, we have also observed that the EPM needs to be used carefully. For configurations of a large number of close-together channels, we discovered during our test runs that the linear system for the time derivatives of the channel parameters (given in appendix B) can become singular, with a rank deficiency of 1. When this occurs during a simulation, the vector of parameters describing the instantaneous configuration can change by an arbitrary amount in a certain direction (the nullspace direction). An example of such a singular configuration is shown in figure 11(a), with an example of how it can suddenly reconfigure in figure 11(b). We do not understand the reason for this singularity, but it is important to emphasise that the singularity is not real – we have confirmed that nothing similar occurs in the full numerical simulations – but, rather, it is an artefact of the approximations made in the EPM when it is used outside of its range of validity.

Figure 11. The EPM singularity for close-together channels. For the hexagonal array of three rings of channels with channel radii 0.054945 and radial separation 0.138756 (a), the linear system for the time derivatives of the channel parameters is singular. As a result, during a simulation, the configuration can ‘jump’ by an arbitrary multiple of the nullspace (b), which is not physically realistic.

This singularity appears to occur in a variety of geometries with a large number of channels whenever the channels are sufficiently close together. It does not appear to be possible to identify in advance whether a given initial configuration will become singular. Because of this, when applying the EPM to geometries with close-together channels, we recommend that the EPM solution be checked to ensure that a singularity has not been crossed. One option is to check the time derivatives for any sudden jumps. Alternatively, the simulation can be run from the final configuration with time reversed and the result compared with the initial configuration. If using a fixed time step, it should be small enough to ensure that the singularity has not been skipped over.

In summary, the EPM agrees well with full numerical simulation within its domain of validity, with identifiable problems occurring outside it. The EPM is therefore a reliable, accurate and fast finite-dimensional reduction of the full infinite-dimensional free-boundary problem in the cross-plane.

Several numerical simulations of fibre drawing have appeared in the literature (Xue et al. Reference Xue, Large, Barton, Tanner, Poladian and Lwin2005a ,Reference Xue, Tanner, Barton, Lwin, Large and Poladian b ,Reference Xue, Tanner, Barton, Lwin, Large and Poladian c ; Chakravarthy & Chiu Reference Chakravarthy and Chiu2009), mainly for MOFs having small numbers of channels. Such approaches are computationally expensive for the complex structures typically seen in MOFs. In contrast, the EPM can be readily solved without the need for effecting a full numerical simulation. This significantly reduces the required computational resources, with the saving becoming more significant, and valuable, as the number of channels increases.

9. Regularisation of the inverse problem

Beyond its advantages in reducing the computational demands for the ‘forward problem’, the generalised EPM has a second, arguably more important, advantage: it provides a stable means of finding solutions of the inverse problem discussed earlier. Simply running the full numerical simulation backwards in time rapidly leads to the growth of small numerical inaccuracies that pollute the computation and lead to unrealistic initial geometries. The EPM model can, on the other hand, be successfully run backwards in time – at least until the geometry becomes unphysical, usually by dint of overlapping channels. In this way, the EPM generates physically realistic initial profiles that are viable candidates for realisation as an actual preform. The EPM therefore provides a natural physically based ‘regularisation’ of this ill-posed inverse problem. (Mathematically, the EPM provides a rational way of remaining close to a ‘slow manifold’ with polluting higher-order modes filtered out.)

In this section, we illustrate this calculation for several example target geometries. In all examples to follow, we use a surface tension ${\it\gamma}=0.25~\text{N}~\text{m}^{-1}$ , typical for silicate glasses (Boyd et al. Reference Boyd, Ebendorff-Heidepriem, Monro and Munch2012). Rather than simply fixing a value of the draw ratio $D$ , we instead fix the outer fibre diameter at the top of the draw to be 1.02 cm and that at the bottom of the draw to be $160~{\rm\mu}\text{m}$ ; these are values typically used in real fibre draws. We fix these diameters in order to more closely emulate what is done in practice: the outer diameter at the top of the draw is determined by the preform diameter while the outer diameter at the bottom of the draw is the desired fibre diameter. The draw ratio $D$ can then be easily calculated from the preform and fibre geometries.

In figures 12(a–d), 13(a–d) and 14(a–f) we pick a target fibre geometry and run the EPM backwards in ${\it\tau}$ time. This allows us to explore possible preform geometries capable of producing this target by picking different values of ${\it\tau}_{L}$ , i.e. the duration of the EPM calculation. The parameter ${\it\tau}_{L}$ should be viewed as ‘labelling’ the possible initial geometries. Because the EPM can be run in either direction, the target fibre geometry is recovered by running any of these initial geometries forwards.

Figure 12. Simple fibre geometry with six channels (a), with three example preform geometries that can be used to produce this fibre: (b)  ${\it\tau}_{L}=0.020$ , (c)  ${\it\tau}_{L}=0.050$ , (d)  ${\it\tau}_{L}=0.075$ . Each preform geometry is obtained by starting with the fibre geometry and running the EPM backwards over a time given by ${\it\tau}_{L}$ . This ${\it\tau}_{L}$ is the reduced time over which sintering will deform that preform into the desired fibre. (e,f) Values of the draw ratio $D$ and physical draw tension ${\it\sigma}$ necessary to draw a preform corresponding to ${\it\tau}_{L}$ into the desired fibre. The values for the examples chosen in the upper plots are marked with crosses.

Figure 13. A more complicated fibre geometry, with preform geometries calculated by running the EPM backwards in the same manner as in figure 12. The same quantities are plotted.

Figure 14. A more realistic fibre geometry (a), with preform geometries calculated by running the EPM backwards, like in figures 12 and 13: (b)  ${\it\tau}_{L}=0.010$ , (c)  ${\it\tau}_{L}=0.020$ , (d)  ${\it\tau}_{L}=0.030$ , (e)  ${\it\tau}_{L}=0.035$ , (f)  ${\it\tau}_{L}=0.040$ . (g,h) Values of the draw ratio $D$ and physical draw tension ${\it\sigma}$ necessary to draw a preform corresponding to ${\it\tau}_{L}$ into the desired fibre. The values for the examples chosen in the upper plots are marked with crosses.

With ${\it\tau}_{L}$ chosen for each example, both the initial and final geometries are known, and so the draw ratio $D$ can be calculated. Equations (3.6b ) and (3.7) then give the value of the tension ${\it\sigma}$ required for this draw. In figures 12(e,f), 13(e,f) and 14(g,h), $D$ and ${\it\sigma}$ are plotted as functions of ${\it\tau}_{L}$ . These are the parameters of interest to experimentalists. Figure 15 shows the profile of the draw for one of the preforms in figure 14.

Figure 15. The profile of the draw for the preform shown in figure 14(f). The configuration at the bottom is the fibre geometry from figure 14. Equation (3.4) is used to obtain the relation between the reduced time and the axial coordinate, with ${\it\sigma}^{\ast }$ obtained by solving (3.6a ).

In these examples, it is evident that the preform configurations corresponding to higher ${\it\tau}_{L}$ are more different from the target fibre configuration – the channels are larger, closer together and more elliptical, characteristics that might make them more difficult to fabricate. On the other hand, the draw tension ${\it\sigma}$ decreases significantly with increasing ${\it\tau}_{L}$ , meaning that these configurations would be easier to draw. This suggests that the best choice of initial configuration will be a compromise between a preform that is easy to fabricate and a tension that is easy to achieve.

How well does the EPM solve the inverse problem? Starting with the configuration from figure 14(f), we compared the forward-time evolution of this preform obtained using the EPM with that given by a full numerical simulation. Figure 16 shows snapshots of both. The full numerical simulation gives a final configuration that closely approximates the final EPM configuration, which is the original fibre geometry from figure 14. Thus, the example preform from figure 14 produces a fibre whose channels are at the desired positions and which closely approximate the target circles.

Figure 16. Snapshots of the forward-time evolution of the preform shown in figure 14(f) obtained using the EPM (solid lines) and using the numerical method (dashed lines): (a)  ${\it\tau}=0$ , (b)  ${\it\tau}=0.020$ , (c)  ${\it\tau}=0.040$ . The middle and right-hand columns show close-ups at each time.

Experimentally, complex preform geometries are becoming easier to achieve with the growing availability of procedures such as the three-dimensional printing of dies (Ebendorff-Heidepriem et al. Reference Ebendorff-Heidepriem, Schuppich, Dowler, Lima-Marques and Monro2014). It is important to remark that the EPM has no theoretical limit on the number of channels – the addition of a new channel simply adds 12 more (real) ordinary differential equations to the system outlined in appendix B.