Industrial dry spinning processes: algorithmic for a two-phase fiber model in airflows

The dry spinning of fibers can be described by three-dimensional multi-phase flow models that contain key effects like solvent evaporation and fiber-air interaction. Since the direct numerical simulation of the three-dimensional models is in general not possible, dimensionally reduced models are deduced. We recently developed a one-two-dimensional fiber model and presented a problem-tailored numerical solution strategy in Wieland et al. (J Comput Phys 384:326–348, 2019). However, in view of industrial setups with multiple fibers spun simultaneously the numerical schemes must be accelerated to achieve feasible simulation times. The bottleneck builds the computation of the radial concentration and temperature profiles as well as their cross-sectionally averaged values. In this paper we address this issue and develop efficient numerical algorithms.


Introduction
In dry spinning processes a polymer solvent mixture is extruded from multiple nozzles into a spinning duct with a tempered air stream. The air stream leads to evaporation of the solvent and leaves solidified fibers behind. At the bottom of the spinning duct the fibers are usually drawn down by a take up roller.
The dry spinning of fibers can be modeled by three-dimensional continuum mechanics. The fibers can be treated as two-phase media of two mutually-penetrating continua (polymer and solvent) interacting with a surrounding aerodynamic flow phase, [13]. However, governing the key effects, like solvent evaporation and interaction between multiple slender fibers and air [21], rises the computational complexity of the multiphase-multiscale problem such that the direct numerical simulation of the three-dimensional models is in general not possible. Therefore, several surrogate fiber models for dry spinning have been developed in the last decades. Starting from uni-axial stationary one-dimensional models based on cross-sectionally averaged balances in [14,15], a two-dimensional model with concentration-dependent rheological laws covering the radial diffusion effects was employed in [4]. A comparative theoretical and experimental study of a cellulose acetate/acetone dry spinning system was presented in [17]. The works of [8,9] extended the models from viscous to viscoelastic material behavior. In [19] we recently developed a dimensionally reduced fiber model which consists of one-dimensional equations for tangential fiber velocity and stress and two-dimensional advection-diffusion equations covering the cross-sectional variations of polymer mass fraction and fiber temperature. We extended this uni-axial fiber model to curved fibers utilizing the theory of Cosserat rods in [20]. Due to the dimensional reduction, certain physical phenomena are certainly not shown by the surrogates. With our model, a description of the shark skin formation is for example not possible due to the neglect of radial velocity information, but the essential effect of solvent evaporation is correctly covered. In comparison to the referential solution of a three-dimensional model, the results are convincing, providing a good approximation, while drastically reducing the computational time, cf. [19]. The numerical solution framework proposed in [19] allows the incorporation of mutual fiber-air interaction in a two-way coupling and, hence, makes the simulation of industrial dry spinning feasible. The simulation results for a cellulose acetate/acetone dry spinning setup stand in accordance to the data in [17]. In view of process design and optimization of complex industrial setups with multiple hundred to more than one thousand fibers spun simultaneously, however, the numerical treatment is still too time-consuming. When solving the coupled one-twodimensional problem iteratively, it turns out that the computation of the radial profiles of polymer mass fraction and fiber temperature as well as their associated cross-sectionally averaged values builds the bottleneck. In this paper, we address this computational challenge. We develop new algorithmic strategies that represent significant improvements on the existing state of the art. Based on a generalized problem description we give details for an efficient implementation of the underlying routines and demonstrate their performance for an industrial example.
The scientific novelty of this paper is the algorithmic; the model and the industrial example are taken from [19]. With regard to a closed presentation and independent readability, the paper is structured as follows. For a better understanding of the underlying multiscalemultiphase problem we start with the three-dimensional fiber description based on a mixture model ansatz in Sect. 2. The dimensionally reduced one-two-dimensional fiber model according to [19] is given in Sect. 3. The core of the paper is Sect. 4, where we introduce new numerical strategies for the computation of the two-dimensional profiles and averaged quantities. Moreover, based on a general problem formulation we describe the coupling concept for the solution of the underlying one-two-dimensional problem and give details for its implementation in view of efficient simulations. In Sect. 5 we apply the framework to the industrial dry spinning of a cellulose acetate-acetone mixture in air. Investigating the computational performance, we discuss the improvements of our developed numerical schemes compared to [19].

Underlying three-dimensional multiphase-multiscale problem
For a better understanding of the dimensionally reduced one-two-dimensional fiber model given in [19], we describe the underlying (Euler-)stationary three-dimensional dry spinning model for a single viscous uni-axial fiber in a surrounding airflow in this section. Since in dry spinning processes no relevant transient effects occur, we restrict to a steady consideration. Let D ⊂ R 3 be the a priori unknown fiber domain whose boundary ∂D = Γ in ∪ Γ fr ∪ Γ out consists of the fixed inlet at the nozzle Γ in , the free lateral fiber surface Γ fr and the outlet Γ out . Proceeding from balance laws for mass, momentum and energy for the two phases i in D, i ∈ {p, d} (polymer and diluent) [13], we employ the mixture model ansatz of [11] to reduce the balances for momentum and energy for the single phases to balances for the mixture. The quantities for the single phases are indicated by the respective index i, i ∈ {p, d}. Note that throughout this paper we use the terms 'diluent' and 'solvent' synonymously.
Let ρ i [kg/m 3 ] and v i [m/s] be the partial densities and velocities for the polymer and diluent phases as well as h i [J/kg] the partial specific enthalpies of polymer and diluent in the mixture, i ∈ {p, d}. Assuming an Eulerian (spatial) description, the stationary mass, momentum and energy balances for the two phases are modeled as  [11] treats the two phases as interpenetrable continua. Its idea is to consider only one momentum equation as sum of the phase balances (1b) and only one energy balance equation as sum of (1c). The mixture density ρ, mixture stress tensor Σ , total body force f, mixture specific enthalpy h as well as the mixture thermal conductivity C are the sums of the quantities of the single phases, i.e., For the mixture we assume ideality, i.e., the volume does not change under mixing and the enthalpy of mixing is zero. This leads to the relations where ρ 0 i , h 0 i denote the material densities and specific enthalpies of pure polymer and solvent. These material parameters might be specified depending on the temperature. Together with the definition of the mixture specific enthalpy h the relation for the enthalpies results in the more common expression The temperature derivatives of h 0 p , h 0 d and h are in particular the specific heat capacities q 0 p , q 0 d and q [J/(kg K)] for constant pressure, yielding Here, we used the fact that the mass fractions of the single phases do not explicitly depend on the temperature, i.e., Moreover, summing up the momentum phase balances (1b) and neglecting diffusive parts in the stresses yields the mixture momentum balance Analogously using Fick's law and the phase balances (1c) we obtain the total energy balance for the mixture Introducing appropriate boundary conditions, the stationary free boundary value problem (BVP) for the fiber unknowns ρ p , p, v, T and D is given by System 1 (Three-dimensional free BVP) Balance laws in D: Kinematic, dynamic, mass and heat flux respective boundary conditions on Γ fr : v · n = 0,

Inlet boundary conditions on
Outlet boundary condition on Γ out : Constitutive laws: Here, we consider body forces due to gravity, i.e., f = ρg, as well as surface forces f [Pa] due to the surrounding airflow. The geometry is specified via the kinematic boundary condition on Γ fr with unit outer normal vector n. At the lateral fiber surface the diluent density has a jump due to the solvent evaporation, moreover, it changes rapidly in the boundary layer that the surrounding air forms around the fiber. The diluent mass flux j c [kg/(m 2 s)] in the aerodynamic boundary layer can be modeled by the difference of the diluent density in the air at the fiber surface ς [kg/m 3 ] and away from the fiber ρ d, [kg/m 3 ] with convective mass transfer coefficient β [m/s], i.e., j c = β(ςρ d, ). Let c = ρ p /ρ = 1ρ d /ρ be the polymer mass fraction, then we particularly use a formulation for j c in terms of the mass fraction associated transfer coefficient γ = β [kg/(m 2 s)] with = ς ρ/ρ d [kg/m 3 ] in System 1. The temperature is continuous at the fiber surface, whereas the heat flux has a jump because of the heat exchange in the air due to the solvent evaporation with evaporation enthalpy δ [J/kg] of the diluent. In the aerodynamic boundary layer the heat flux j T [W/m 2 ] is described-analogously to the mass flux-by the difference of the temperature at the fiber surface and away from the fiber T [K] with heat transfer coefficient α [W/(m 2 K)]. The parameters δ and ς might be functions of c and T, whereas the transfer coefficients α and β depend on the state of the surrounding airflow and especially on the relative velocity between fiber and airflow. In System 1 the surrounding airflow is assumed to be known in the sense that the quantities f , ρ d, and T as well as the airflow dependencies of α and β are given for each point of the fiber surface. Moreover, the parameters D, C and μ might depend on the mass fraction c and temperature T with suitable models assumed to be available.
Remark 1 (Mass transfer) The used formulation of the diluent mass transfer j c is motivated from the fact that at the fiber surface the diluent density in the air is mainly linearly proportional to the diluent density in the fiber, i.e., ς ∼ ρ d . Setting ς = (ρ d /ρ) separates the linear part from the remainder . The mass transfer is principally driven from the linear part in terms of a Robin-type boundary condition, whereas the remainder is incorporated in the transfer coefficient γ = β . This splitting might allow a different treatment of the terms, which becomes essentially for our numerical treatment of the problem.

Reduced fiber dry spinning model
Due to the complexity of the three-dimensional problem (System 1) that prevents the simulations of industrial dry spinning setups, we proposed a dimensionally reduced model in [19] that is of good approximation quality and efficiently evaluable. It combines the relevant two-dimensional aspects (radial profiles for polymer mass fraction c and temperature T), that are essential due to the evaporation of solvent [14,15], with one-dimensional cross-sectionally averaged balances for tangential fiber velocity u and stress σ . In dimensionless form the model for u, σ , c, T is given by System 2 (One-two-dimensional BVP) Cross-sectionally averaged balance laws, s ∈ (0, 1): with boundary conditions at inlet s = 0 and outlet s = 1: Radial equations, (r, s) ∈ (0, 1) 2 : with boundary conditions at inlet s = 0, fiber surface r = 1 and symmetry boundary r = 0: Constitutive laws and geometric relation: Abbreviations: The one-dimensional equations for the tangential fiber velocity u and stress σ result from averaging the Newtonian stress tensor Σ and the momentum balance in System 1 over circular fiber cross-sections. The two-dimensional advection-diffusion equations for the polymer mass fraction c and temperature T, which reveal the radial effects that are essential due to evaporation, are obtained from the diluent mass and energy balance in System 1 by the assumption of radial symmetry and by linearization around the crosssectional averaged polymer mass fractionc and temperatureT. Moreover, averaging the polymer mass balance in System 1 over the fiber cross-sections yields the constant polymer mass flux Q =cρ(c,T)uR 2 π = c in ρ(c in , T in )u in R 2 in π . This means we do not face a free boundary value problem anymore, but can conclude the fiber radius R from the boundary conditions at the inlet. For further details see [19].
Remark 2 (Non-dimensionalization and scaling) System 2 is formulated with respect to dimensionless quantities y that have been scaled onto the space domain (0, 1) (onedimensional quantities) and (0, 1) 2 (two-dimensional quantities), respectively. The nondimensionalization and scaling are done in separate steps. First, dimensionless quantities are introduced as y(s) =y(s 0s )/y 0 ,ỹ(r,s) =y(d 0r , s 0s )/y 0 for any scalar-or vector-valued dimensionful quantityy and corresponding reference value y 0 . Due to technical reasons we use slightly different scales compared to [19], for the reference values and resulting dimensionless numbers see Table 1. After nondimensionalization the equations are considered on the non-dimensional space domain Table 1 Overview over composite reference values used for non-dimensionalization of the uni-axial fiber dry spinning model and the resulting dimensionless numbers. Here, the scales M,0 , v 0 , r 0 , d 0 ,

Dimensionless numbers
Description Formula (0, L) (one-dimensional equations) and (0, R(s)) × (0, L) (two-dimensional equations), respectively. Second, to simplify the numerical treatment of the equations, one-dimensional equations are transformed onto the space domain (0, 1) and two-dimensional equations onto the space domain (0, 1) 2 via respectively. Since the fiber radius R is an a priori unknown, this domain transformation omits domain changes during the numerical solution of the problem, which would cause computational expensive re-meshing of the computational domain. The transformation inserts the unknown R into our model equations, but this can efficiently be handled with our numerical solution algorithm. Since the dimensionful fiber lengthL is a priori known, choosing s 0 = r 0 =L would lead to a fiber length equal to one in non-dimensional form (i.e., L = 1).
Remark 3 (Averaged viscosity) In System 2 we could approximate the cross-sectionally averaged dynamic viscosity in the same way as the densities, i.e., However, due to in general highly nonlinear rheological models for the dynamic viscosity, this approximation is not feasible for industrial setups.
As aerodynamic force model f air we use with normalized fiber tangent τ , τ = 1, and dimensionless drag function F : S 2 × R 3 → R 3 given in [12], where S 2 denotes the three-dimensional sphere. The relative velocity between fiber and airflow reads v rel = vuτ . Note that to distinguish the fiber quantities from the airflow quantities all airflow associated fields are labeled with the index . In particular, v denotes the velocity, ρ the density, ν the kinematic viscosity, λ the thermal conductivity, and q the specific heat capacity of the air. Moreover, D d, denotes the diffusivity of diluent in the air and ρ d, the diluent density in the air away from the fiber. All these quantities are space-and time-dependent fields assumed to be dimensionless and known-for example provided by an external computation. The corresponding reference values used for non-dimensionalization are denoted with the index 0 and given in Table 1. Furthermore, we employ the models for the heat and mass transfer with the associated dimensionless function N : R 3 → R and the model for the diluent density in the air at the fiber surface given in [19]. The closing of System 2 requires further models for the diffusion coefficient D, thermal conductivity C, evaporation enthalpy δ, specific heat capacities q 0 p , q 0 d , dynamic viscosity μ as well as material densities ρ 0 p , ρ 0 d . For our industrial setup in Sect. 5 we employ detailed models taking the fiber materials and rheological effects into account.

Numerical framework
The fiber model for dry spinning of a single straight fiber (System 2) is a coupled system of one-and two-dimensional model equations. While the one-dimensional equations (2a) can be formulated as parametric boundary value problem of ordinary differential equations, the two-dimensional equations (2b) form a system of radial advection-diffusion equations with Robin-type boundary conditions. In [19] we presented a numerical framework for these two subproblems together with an iterative coupling procedure that made simulations of dry spinning processes feasible. However, when multiple hundred fibers are spun simultaneously and a two-way coupling with the surrounding airflow has to be taken into account, the simulations are still very time-consuming and can easily involve unfeasible computation times for industrial setups. When solving the coupled one-twodimensional problem, it turns out that the computation of the two-dimensional profiles of polymer mass fraction and fiber temperature forms the bottleneck. In this section we, therefore, develop an algorithmic procedure to speed up the computation of these quantities. The key idea is that no radial profiles have to be computed in scenarios, where only surface and averaged values of the two-dimensional quantities are needed. The associated computation then simplifies to the solution of an integral equation and an ordinary differential equation. In the numerical treatment the solution of the integral equation reduces to the solution of a linear system of equations with lower triangular system matrix, such that it can be solved by forward substitution. When two-dimensional profiles are needed, e.g., for the dynamic viscosity (cf. Remark 3), we propose a new algorithm. In the following we give a general description of the coupled problem, present our new numerical strategies, describe the strategy for the iterative solving of the coupled problem, and give a survey for the further weak iterative coupling with an airflow simulation.

General problem formulation
To illustrate the structure of System 2, we embed it into the class of problems that consist of a one-dimensional boundary value problem of the form The diffusion coefficient and boundary conditions in the radial advection-diffusion equation in (3a) depend on the variable ψ such that the associated problem is nonlinear. This nonlinearity is tackled by decoupling into the following two problems for given functions y,ψ . In particular, we note that (3c) is a linear problem. Subsequently, we explain the respective numerical algorithms for the decoupled problems (3b), (3c) and discuss a coupling strategy to obtain a solution of the coupled problem (3a).
In view of the coupling we need the associated radial advection-diffusion equations averaged over the fiber cross-sections, which we already introduce here
Here, 1 ∈ R n denotes the n-dimensional tuple of ones. The functions f 0 , g 0 are chosen in such a way that for p = 0 an analytical solution is known. Given this starting solution, we seek for a sequence of parameter tuples 0 = p 0 , p 1 , . . . , p m = 1 such that the solution to the respective predecessor boundary value problem provides a good initial guess for the successor. The solution associated to p = 1 finally belongs to the original system. With the help of the continuation parameters certain terms in the ordinary differential equation can be first excluded, then included. Also the boundary conditions can be varied. The core of a robust continuation procedure are the choice of a continuation path and the step size control to navigate through a high-dimensional parameter space. They decide about failure or success because there are not always existing solutions and several ways might be possible. For the choice of the continuation path we refer to [19] and the step size control follows the strategy given in [2].

Radial advection-diffusion equation algorithm
We consider (3c) with given functions y,ψ . This means λ(y(s),ψ(·, s), s) = λ(s), such that the problem becomes linear. Furthermore, we transform the unknown ψ by the following rule  , x), as well as functions a, b : [0, 1] → R and constant φ in ∈ R. If the function a is constant, the solution can be given analytically in terms of an explicit expression. For non-constant a we treat the Robin boundary condition as Neumann boundary condition with right hand side depending on φ and find an implicit solution expression with Green's function g in [3,7,16]. The solution expression reads where J i denote the i-th Bessel functions of first kind and β m > 0, m ∈ N, are the (nontrivial) ascending zeros of the first Bessel function of first kind, i.e., J 1 (β m ) = 0. These values are tabulated in literature, see e.g. [7]. For φ| r=1 the solution expression yields a Volterra integral equation of second kind The integral kernel g is singular for x = 0. Therefore, numerical integration in the sense of quadrature formulas cannot be applied directly to the integral in (7), because they involve the evaluation of the integrand function at or close to the singularity. Hence, we use the product integration method that we introduced in [18,19]. This means we substitute the function k (·, φ(1, ·)) piecewise by Lagrange interpolation polynomials and employ iterated integration by parts to isolate the singularity of the kernel function g. Here, we choose constant polynomials corresponding to the implicit Euler method yielding an A-and Lstable method. For the discussion of stability properties as well as an extension to a higher order method with quadratic polynomials and nodes corresponding to the Lobatto IIIa collocation scheme we refer to [19]. Let 0 = x 0 < x 1 < · · · < x N x = 1 be the mesh points in x-direction. Substituting the integrand function k (·, φ(1, ·)) piecewise by constants in the sense of the implicit Euler scheme and subsequent integration yields which we evaluated analytically. This results in the linear system of equations for the unknown vector φ(1) = (φ 1 (1), . . . , φ N x (1)) ∈ R N x . Here, φ in = φ in 1 ∈ R N x with 1 ∈ R N x the N x -dimensional tuple of ones and I ∈ R N x ×N x denotes the identity matrix. The matrix G(r) = (G ij (r)) i,j=1,...,N x ∈ R N x ×N x is defined as else.
Furthermore, we use the vectors a = (a 1 , . . . , The system matrix (I -2πG(1) · diag(a)) is lower triangular, such that the linear system of equations can be solved by forward substitution. After solution of the integral equation for the boundary values φ(1) the profile φ(r) = (φ 1 (r), . . . , φ N x (r)) ∈ R N x , 0 ≤ r ≤ 1, can be calculated at the mesh points x i using the same approximation strategy for the singular integral For the calculation of the matrix G(r), we define further lower triangular matrices H(β m ) = (H ij (β m )) i,j=1,...,N x ∈ R N x ×N x and = ( ij ) i,j=1,...,N x ∈ R N x ×N x as Thus, we can re-write the entries of G(r) , and the profiles are given as Due to the separation of the r-dependency and its exclusive occurrence in the factor B, this representation allows an efficient implementation of the calculation of the profiles with complexity in O(max{N 2 x N β , N r }) with N β the number of zeros β m of J 1 taken into account for the calculation. This is a great improvement compared to the implementation given in [19] that has the complexity O(N 2 x N β N r ). In view of System 2 the cross-sectionally averaged quantityφ is needed. This quantity can be obtained with the help of two different approaches • Method 1: Averaging of the profiles, • Method 2: Solution of the integrated advection-diffusion equation.
In Method 1 the calculation follows directly from the profiles Introducing a mesh in r-direction, 0 = r 0 < · · · < r N r = 1, N r ∈ N, and utilizing the trapezoidal rule we get The complexity of Method 1 is O (N x N r ).
In Method 2 the averaged quantityφ is calculated with the help of the integrated advection-diffusion equation (4) that reads component-wisely after transformation onto the x-grid This differential equation can be solved by integration with respect to x omitting any discretization error in radial direction. Employing right Riemann summation, which corresponds to the product integration method based on the implicit Euler scheme, we get the recursion formulā In difference to Method 1, Method 2 omits any radial resolution and has the lower complexity of O(N x ). Moreover, in situations where only boundary φ| r=1 and averaged valuesφ are needed the calculation of two-dimensional profiles φ can be omitted completely when Method 2 is used. Thus, Method 2 is obviously preferable to Method 1 due to efficiency reasons.

Coupling procedure
For the solution of (3a) we proposed an iterative coupling of the solution algorithm for the one-dimensional boundary value problem with the solution algorithm for the radial advection-diffusion equation in [19]. Here, we formulate this coupling procedure with respect to the general problem formulation given in Sect. 4.1.
In view of the radial advection-diffusion equation a solution of the integrated problem (4) based on averaged quantities serves as initial guess in our coupling procedure. Since we employ the profiles of polymer mass fraction c and temperature T to approximate the averaged dynamic viscosity μ(c, T) R 2 /(πR 2 ) (cf. Remark 3) we have to build up the two-dimensional profiles of these quantities. In total the coupling of the one-and twodimensional model equations is done as described in Algorithm 1.
In the coupling procedure (Algorithm 1) cross-sectionally averaged quantities are needed for the solution of the one-dimensional boundary value problem in Line 5, i.e., cross-sectionally averaged quantities propagate into the operator B. These averaged quantities have to be calculated after the computation of the two-dimensional quantities ψ i+1 in Line 4. This can be done by either of the Methods 1 or 2. In view of its implementation Method 2, which involves the solution of the integrated equation (4), needs further considerations. We introduce two submethods for the computation of the averaged quantitȳ ψ i+1 : This means, Method 2a utilizes the profile information from the preceding iteration step, whereas Method 2b uses the updated information from the actual iteration step. At first glance, Method 2b looks more straightforward but it has a major drawback: The nonlinearity of the coupled problem (3a) comes from the nonlinearity of the one-dimensional and two-dimensional equations themselves as well as the nonlinear coupling of these equations. Whereas the nonlinearity of the one-dimensional equations is treated by the continuation method, the remaining nonlinearities in common are tackled iteratively by Algorithm 1. This means, ψ i+1 as calculated in Line 4 is only a solution of ψ i+1 -D[y i , ψ i+1 ] when Algorithm 1 is converged. Therefore, when solving the integrated equation (4),ψ Remark 4 (Two-way fiber-air interaction) Since the nature of the airflow in the spinning duct has a deep impact onto the behavior of the fibers and also the fibers itself affect the properties of the airflow, a two-way coupling between the fibers and the airflow has to be taken into account. This coupling is realized by homogenized exchange models. The exchange models are constructed in such a way, that the principle of action equaling reaction is fulfilled, cf. [5,6]. The two-way fiber-air-interaction is realized by a weak coupling algorithm which iterates between fiber computations and airflow simulations. This procedure allows the combination of self-implemented code for the fiber dynamics as well as a commercial software for the air dynamics and was successfully used in studies on glass wool manufacturing [2]. In total, our complete procedure consists of nested iterations: An inner iteration for the coupling of one-and two-dimensional fiber equations (Algorithm 1) and an outer iteration realizing the mutual fiber-airflow interaction. For further details we refer to [19].

Industrial application
In this section, we consider the industrial example of [19] and investigate the performance of the new numerical methods. In particular, we highlight the efficiency of the calculation of two-dimensional profiles (c, T) using the representation (8) and compare Method 1, Method 2a and Method 2b for the computation of cross-sectionally averaged quantities (c,T ). Note that for the sake of clarity we use dimensionful quantities to describe the following setup and results.
Remark 5 Algorithm 1 is implemented in MATLAB a version R2016b, where the BVP solver bvp4c.m is used with the default values. For the computation of the twodimensional profiles we choose N r = 201 equidistant points in r-direction and N β = 318, i.e., we use all non-trivial zeros of J 1 smaller than 10 3 . As stopping criterion for the iteration in Algorithm 1 we use tol = 10 -6 . For the simulation of the airflow we employ ANSYS Fluent b with a solver accuracy of O(10 -6 ). The break-up criterion of the fiber-air coupling algorithm satisfies an error tolerance tol = 10 -5 . For the considered industrial setup the fiber-air coupling algorithm converges after 30 iteration steps.
In the setup 60 fibers are spun in a spinning duct and drawn down by a take up roller at the duct bottom. The length L and tangent vector τ of a fiber are given through its inflow position in the spinneret and the position of the take up at the bottom. We consider the outer basis {a 1 , a 2 , a 3 } ⊂ R 3 as given in Fig. 1. The holes in the spinneret are ordered in three parallel rows with 20 holes each. In particular, we consider the fibers of the middle row to be spun in the x 1 -x 3 -plane corresponding to x 2 = 0. The airflow forms a counter-current to the fiber flow and we assume dry air at the inlet (pipe radius R ,in = 9.8· 10 -2 m) with an absolute air inflow velocity (in negative a 1 -direction) at the bottom pipe of v = 0.22 m/s and air inflow temperature T = 330 K. A cellulose acetate (CA)-acetone mixture is dry spun. For the closing of our dry spinning model (System 2) we employ rheological models for the quantities C, D, h 0 d , δ, μ that are given in [19]. All further process parameters, the specific reference values as well as the dimensionless numbers are listed in the Appendix.
Cellulose acetate/acetone dry spinning was previously also investigated in [9,17]. The results obtained by our model-simulation framework in [19] stand in accordance to the data in [17]. Exchanging the algorithmic by the new numerical methods (developed in Sect. 4) has no effect on the approximation quality, but only on the computational effort of the simulation. However, before we discuss the performance of the new methods, we present the simulation outcome to illustrate the typical fiber behavior for an industrial dry spinning setup. The solutions of the 60 fibers do not show visible differences, such that we illustrate-as an example-the solution behavior for the fiber that is centrally located in the spinning duct (cf. Fig. 1) in Fig. 2. The cross-sectionally averaged CA mass fraction increases from c 0 = 0.29 at the inlet toc = 0.90 at the outlet, indicating the evaporation of the acetone during the spinning process. The fiber speed u reaches the take up speed u out = 10 m/s approximately 1 m away from the nozzle. The tensile stress σ increases correspondingly to the changes of u and the integrated viscosity μ R 2 . Directly at the nozzle the fiber loses heat due to acetone evaporation and away from the nozzle the fiber temperature approaches the airflow temperature. The fiber thinning is caused by the fiber take up as well as the evaporation of diluent. The final fiber radius is R = 1.06 · 10 -5 m (cf. Fig. 3). The effect of solvent diffusion in the fiber interior can clearly be seen at the profile for the polymer mass fraction in Fig. 3 showing an inhomogeneous CA-acetone distribution in radial direction. Since there are no visible radial effects in the temperature due to the fact that the inverse of the temperature associated Peclet number 1/Pe T is three magnitudes larger compared to the mass associated one 1/Pe c (cf. Table 4) we omit to show a radial temperature profile. Investigating the performance of the new numerical methods, first, we highlight the efficiency of the calculation of profiles using the representation (8). Whereas in [19] the total simulation time with a straightforward implementation for the profile solution (regardless of (8)) was 7.95 h, the simulation only needs 5.20 h utilizing the efficient profile calculation (cf. (8)) here, which means a reduction of 34.6% in total. Exemplarily considering N r = 201 and N x = 501 points in r-and x-direction for the profiles as well as N β = 318, the calcu- Second, we compare the three different procedures (Method 1, Method 2a, and Method 2b) for the computation of the cross-sectionally averaged polymer mass fraction c and consider Method 1 (averaging of two-dimensional profiles), which has been used in [19], as reference. We exemplarily consider the last step in the air-fiber coupling (30th iteration) where Algorithm 1 converges after 13 iteration steps. Figure 4 shows the results after the 2nd and the 13th iteration step. We clearly see that Method 2b initially produces an inconsistent solution due to the use of inconsistent profile information, see Sect. 4.4. The calculated averaged polymer mass fractionc exceeds the average of the associated profile (solution obtained by Method 1) greatly. This behavior could lead to a completely unphysical solution behavior in the next iteration steps and even to a failure of the iterative coupling procedure (divergence of Algorithm 1). On the other hand, Method 2a yields a consistent solution. As soon as Algorithm 1 is converged, the solutions obtained with Method 1, Method 2a and Method 2b coincide. This fact can be seen in Fig. 4 indicating the relative L 2 -error between these three solutions with Method 1 as reference. Consequently, Method 2a should be used for the calculation ofc,T. In view of computational time the difference between Method 1 and Method 2a is marginal. Whereas Method 1 takes 1.30 · 10 -3 s, Method 2a needs 3.42 · 10 -4 s for N x = 501 points. However, the clear advantage of Method 2a comes into play when no two-dimensional profiles need to be calculated.

Conclusion
For the modeling of viscous fibers in dry spinning processes one-dimensional equations for tangential fiber velocity and stress coupled with two-dimensional radial advectiondiffusion equations for polymer mass fraction and temperature are used. In this paper we introduced new algorithmic concepts for the fast computation of the two-dimensional information. The efficient calculation of surface values and averaged quantities as well as the speed up of the coupling procedure for one-and two-dimensional equations yields feasible computation times for simulations where multiple hundred fibers are spun simultaneously and mutual air-fiber interaction effects are incorporated. Consequently, our proposed model-simulation framework provides the possibility for optimizing industrial dry spinning devices. Depending on the mixture to be spun and the device the rheological models as well as the material and process parameters must be certainly adapted.

Appendix: Data to industrial process
In this paper (cf. Sect. 5) we consider the same dry spinning setup that was studied in [19] and demonstrate the performance of the new numerical methods (cf. Sect. 4). We assume that a polymer solvent mixture consisting of cellulose acetate (CA)-acetone is dry spun. Table 2 Physical parameters of the industrial dry spinning setups taken from [9,17] and setup-specific reference values Table 4 Dimensionless numbers for one representative fiber (cf. Fig. 1) in the industrial dry spinning setup assuming the airflow without fibers (first step in the fiber-airflow coupling). Note that some dimensionless numbers might vary during the simulation due to the fiber-airflow interaction (cf. A similar setup has been investigated in [9,17], from where we take the physical material parameters for the polymer CA and the diluent acetone, see Table 2. The setup-specific reference values used for the non-dimensional form of the model equations (cf. System 2) are also given in Table 2. All further process parameters can be found in Table 3. Utilizing the setup-specific reference values (cf. Table 2) the resulting dimensionless numbers are given in Table 4.