Melt-Blowing of Viscoelastic Jets in Turbulent Airflows: Stochastic Modeling and Simulation

TL;DR

This paper develops a stochastic modeling framework to simulate melt-blowing viscoelastic jets in turbulent airflow, revealing turbulence's crucial role in fiber thinning and producing realistic fiber diameter predictions.

Contribution

It introduces a novel stochastic aerodynamic force model incorporating turbulence effects into viscoelastic jet simulations, improving accuracy over previous methods.

Findings

Turbulence significantly influences jet thinning.

Simulated fiber diameters match experimental observations.

The model enables feasible industrial-scale simulations.

Abstract

In melt-blowing processes mico- and nanofibers are produced by the extrusion of polymeric jets into a directed, turbulent high-speed airflow. Up to now the physical mechanism for the drastic jet thinning is not fully understood, since in the existing literature the numerically computed/predicted fiber thickness differs several orders of magnitude from those experimentally measured. Recent works suggest that this discrepancy might arise from the neglect of the turbulent aerodynamic fluctuations in the simulations. In this paper we confirm this suggestion numerically. Due to the complexity of the process direct numerical simulations of the multiscale-multiphase problem are not possible. Hence, we develop a numerical framework for a growing fiber in turbulent air that makes the simulation of industrial setups feasible. For this purpose we employ an asymptotic viscoelastic model for the fiber. The turbulent effects are taken into account by a stochastic aerodynamic force model where the underlying velocity fluctuations are reconstructed from a $k$-$\epsilon$ turbulence description of the airflow. Our numerical results show the significance of the turbulence on the jet thinning and give fiber diameters of realistic order of magnitude.

Figures

Tables

Table 1. Table 2.1. Overview over reference values used for non-dimensionalization and the resulting dimensionless numbers.

Reference values
Description	Formula		Unit
fiber curve	$r_0$ =	$H$	m
fiber diameter	$d_0$ =	$d_{in}\sqrt{\pi }/2$	m
fiber velocity	$v_0$ =	$v_{in}$	m/s
fiber temperature	$T_0$ =	$T_{in}$	K
fiber mass line density	${\varrho }_{M0}$ =	$\rho d_0^2$	kg/m
fiber volume line density	${\varrho }_{V0}$ =	$d_0^2$	m2
fiber stress	${\sigma }_0$ =	${\varrho }_{M0}{v}_0^2/d_0^2$	Pa
fiber pressure	$p_0$ =	${\sigma }_0$	Pa
fiber kinematic viscosity	${\mu }_0$ =	$\mu \left(T_0\right)$	Pas
fiber relaxation time	${\theta }_0$ =	$\theta \left(T_0\right)$	s
outer forces	$f_0$ =	${\varrho }_{M0}{v}_0^2/r_0$	N/m
heat transfer coefficient	${\alpha }_0$ =	${\alpha }_{in}$	W/(m2K)
length scale	${\zeta }_0$ =	$r_0$	m
time scale	$t_0$ =	$r_0/v_0$	s
air velocity	$v_{\star ,0}$ =	$v_0$	m/s
air density	${\rho }_{\star ,0}$ =	${\rho }_{\star ,in}$	kg/m3
air kinematic viscosity	${\nu }_{\star ,0}$ =	${\nu }_{\star ,in}$	m2/s
air specific heat capacity	$c_{p,\star ,0}$ =	$c_{p,\star ,in}$	J/(kgK)
air thermal conductivity	${\lambda }_{\star ,0}$ =	${\lambda }_{\star ,in}$	W/(mK)
air turbulent kinetic energy	$k_{\star ,0}$ =	$k_{\star ,in}$	m2/s2
air viscous dissipation rate	${ϵ}_{\star ,0}$ =	${ϵ}_{\star ,in}$	m2/s3

Table 2. Table 4.1. Overview over process and physical parameters in the industrial melt-blowing setup according to [ 16 ] and the resulting dimensionless numbers.

Parameters
Description	Symbol	Value	Unit
device height	$H$	$1.214\cdot {10}^{-1}$	m
nozzle diameter	$d_{in}$	$4\cdot {10}^{-4}$	m
fiber speed at nozzle	$v_{in}$	$1\cdot {10}^{-2}$	m/s
fiber temperature at nozzle	$T_{in}$	$5.532\cdot {10}^2$	K
heat transfer at nozzle	${\alpha }_{in}$	$1.595\cdot {10}^3$	W/(m2K)
polar angle at nozzle	${\vartheta }_{in}$	$\pi /2$	–
azimuth angle at nozzle	${\phi }_{in}$	$\pi$	–
fiber density	$\rho$	$7\cdot {10}^2$	kg/m3
fiber specific heat capacity	$c_p$	$2.1\cdot {10}^3$	J/(kgK)
end time	$t_{end}$	$2.20\cdot {10}^{-2}$	s
air density at nozzle	${\rho }_{\star ,in}$	$1.187$	kg/m3
air kinematic viscosity at nozzle	${\nu }_{\star ,in}$	$1.8\cdot {10}^{-5}$	m2/s
air specific heat capacity at nozzle	$c_{p,\star ,in}$	$1.006\cdot {10}^3$	J/(kgK)
air thermal conductivity at nozzle	${\lambda }_{\star ,in}$	$2.42\cdot {10}^{-2}$	W/(mK)
air turbulent kinetic energy at nozzle	$k_{\star ,in}$	$2.181\cdot {10}^2$	m2/s2
air viscous dissipation rate at nozzle	${ϵ}_{\star ,in}$	$1.808\cdot {10}^7$	m2/s3

Table 3. Table 4.2. Dimensionless numbers characterizing the fiber behavior in the stochastic region.

Dimensionless numbers
Description	Symbol	Value
slenderness	$\epsilon$	$4.47\cdot {10}^{-5}$
Reynolds	$\mathrm{Re}$	$2.15\cdot {10}^2$
Deborah	$\mathrm{De}$	$2.72\cdot {10}^2$
Froude	$\mathrm{Fr}$	$4.84\cdot {10}^1$
Stanton	$\mathrm{St}$	$1.75\cdot {10}^{-4}$
Mach	$\mathrm{Ma}$	$2.42\cdot {10}^2$
air drag associated	$A_{\star }$	$2.66\cdot {10}^1$
mixed (air-fiber) Reynolds	${\mathrm{Re}}_{\star }$	$1.40\cdot {10}^1$
Nusselt	${\mathrm{Nu}}_{\star }$	$2.68$
Prandtl	${\mathrm{Pr}}_{\star }$	$6.23\cdot {10}^{-1}$
turbulence degree	${\mathrm{Tu}}_{\star }$	$7.34\cdot {10}^{-1}$
turbulent time	${\mathrm{Tt}}_{\star }$	$1.14\cdot {10}^2$

Equations

$${\partial }_t\mathbf{r}$$

$${\partial }_{\zeta }\mathbf{r}$$

$${\partial }_t({\varrho }_M\mathbf{v})$$

$$c_p{\partial }_t({\varrho }_{M}T)$$

$${\partial }_t\sigma$$

$${\partial }_{t}p$$

$$d=2\sqrt{\frac{{\varrho }_V}{\pi \Vert \mathbf{\tau }\Vert }}.$$

$${\partial }_t\mathbf{r}$$

$${\partial }_{\zeta }\mathbf{r}$$

$${\partial }_t\mathbf{v}$$

$${\partial }_{t}T$$

$$\mathrm{De}\left({\partial }_t\sigma -(2\sigma +3p)\frac{{\partial }_t\Vert \mathbf{\tau }\Vert }{\Vert \mathbf{\tau }\Vert }\right)+\frac{\sigma }{\theta }$$

$$\mathrm{De}\left({\partial }_{t}p+p\frac{{\partial }_t\Vert \mathbf{\tau }\Vert }{\Vert \mathbf{\tau }\Vert }\right)+\frac{p}{\theta }$$

$${\partial }_t\mathbf{\phi }+\mathbf{M}(\mathbf{\phi })\cdot {\partial }_{\zeta }\mathbf{\phi }+\mathbf{m}(\mathbf{\phi })=0$$

$$\sigma (0,t)=0,p(0,t)=0.$$

$$e(0,t)=1,$$

$$\mathbf{r}(-{\zeta }_L(t),t)$$

$$\sigma (-{\zeta }_L(0),0)={\sigma }_{in}/{\sigma }_0,p(-{\zeta }_L(0),0)=p_{in}/p_0,e(-{\zeta }_L(0),0)=1.$$

$$\sigma (-{\zeta }_L(t),t)$$

$$e(-{\zeta }_L(t),t)$$

$${\partial }_{t}e-\mathbf{t}\cdot {\partial }_{\zeta }\mathbf{v}=0.$$

$$\mathbf{t}(\vartheta ,\phi )=(\sin \vartheta \cos \phi ,\sin \vartheta \sin \phi ,\cos \vartheta ),\vartheta \in [0,\pi ],\phi \in [0,2\pi ).$$

$${\partial }_t\mathbf{t}={\partial }_t\left(\frac{\mathbf{\tau }}{e}\right)=\frac{1}{e}\left(\mathbf{I}-\mathbf{t}\otimes \mathbf{t}\right)\cdot {\partial }_{\zeta }\mathbf{v}$$

$${\partial }_t\vartheta =\frac{1}{e}\mathbf{n}\cdot {\partial }_{\zeta }\mathbf{v},{\sin }^2\vartheta {\partial }_t\phi =\frac{1}{e}\mathbf{b}\cdot {\partial }_{\zeta }\mathbf{v}.$$

$${\partial }_t\mathbf{r}-\mathbf{v}$$

$${\partial }_{t}e-\mathbf{t}\cdot {\partial }_{\zeta }\mathbf{v}$$

$${\partial }_t\vartheta -\frac{1}{e}\mathbf{n}\cdot {\partial }_{\zeta }\mathbf{v}$$

$${\sin }^2\vartheta {\partial }_t\phi -\frac{1}{e}\mathbf{b}\cdot {\partial }_{\zeta }\mathbf{v}$$

$${\partial }_t\mathbf{v}-{\partial }_{\zeta }\left(\sigma \frac{\mathbf{t}}{e}\right)-\frac{1}{{\mathrm{Fr}}^2}{\mathbf{e}}_g-{\mathbf{f}}_{air}$$

$${\partial }_{t}T+\frac{\mathrm{St}}{\epsilon }\pi d\alpha (T-T_{\star })e$$

$$\mathrm{De}{\partial }_t\sigma +\left(-\mathrm{De}(2\sigma +3p)-\frac{\mu }{\theta }\frac{3}{\mathrm{Re}}\right)\frac{\mathbf{t}}{e}\cdot {\partial }_{\zeta }\mathbf{v}+\frac{\sigma }{\theta }$$

$$\mathrm{De}{\partial }_{t}p+\left(\mathrm{De}p+\frac{\mu }{\theta }\frac{1}{\mathrm{Re}}\right)\frac{\mathbf{t}}{e}\cdot {\partial }_{\zeta }\mathbf{v}+\frac{p}{\theta }$$

$$\mathbf{r}(-{\zeta }_L(t),t)$$

$$\mathbf{v}(-{\zeta }_L(t),t)$$

$$e(-{\zeta }_L(0),0)=1,\sigma (-{\zeta }_L(0),0)={\sigma }_{in}/{\sigma }_0,p(-{\zeta }_L(0),0)=p_{in}/p_0,$$

$$\text{if }{\lambda }_3>-1\text{: }$$

$$\text{if }{\lambda }_5>-1\text{: }$$

$$e(0,t)=1,\sigma (0,t)=0,p(0,t)=0.$$

$$\displaystyle\mathbf{f}_{air}=e\frac{\mathrm{A}_{\star}}{\mathrm{Re}_{\star}^{2}}\frac{\rho_{\star}\nu^{2}_{\star}}{d}\mathbf{F}\bigg{(}\mathbf{t},\mathrm{Re}_{\star}\frac{d}{\nu_{\star}}(\mathbf{v}_{\star}-\mathbf{v})\bigg{)},\qquad\mathbf{F}(\mathbf{t},\mathbf{w})=r_{n}(w_{n})\mathbf{w_{n}}+r_{t}(w_{n})\mathbf{w_{t}}.$$

$$\alpha =\frac{1}{{\mathrm{Nu}}_{\star }}\frac{{\lambda }_{\star }}{d}\mathcal{N}\left({\mathrm{Re}}_{\star }\frac{d}{{\nu }_{\star }}({\mathbf{v}}_{\star }-\mathbf{v})\cdot \mathbf{t},{\mathrm{Re}}_{\star }\frac{d}{{\nu }_{\star }}\Vert {\mathbf{v}}_{\star }-\mathbf{v}\Vert ,{\mathrm{Pr}}_{\star }\frac{c_{p,\star }{\rho }_{\star }{\nu }_{\star }}{{\lambda }_{\star }}\right).$$

Full text

Keywords. melt-blowing, fiber dynamics, upper-convected Maxwell model, turbulence modeling, boundary value problem, finite volume scheme

AMS-Classification. 35Lxx, 68U20, 76-XX

1. Introduction

Melt-blowing is a widely used production method for polymer micro- and nanofibers economically attractive due to low production costs. Fabrics of meltblown fibers are nonwovens, e.g., filters, hygiene products, battery separators. Details on the technology can be found in [10, 24]. A typical setup of a melt-blowing device is illustrated in Fig. 1.1. In the process, molten polymer is fed through a nozzle into a forwarding high-speed and highly turbulent air stream to be stretched and cooled down. The resulting fibers are laid down onto some collector, e.g., conveyor belt. In contrast to melt-spinning processes, where the stretching is caused by a mechanical take-up, in melt-blowing the fiber jet thinning is due to the driving high-velocity air stream with its turbulent nature.

To deepen the understanding on the mechanism of jet thinning in melt-blowing extensive diverse studies have been performed in the last years, covering experimental investigations, e.g., [4, 6, 11, 32, 33], combined experimental numerical works, e.g., [26, 30, 34], as well as numerical computations, e.g., [7, 25, 29, 36]. However, so far, there is an obvious gap between the experimental and numerical results for the achieved fiber thickness in literature. The existing numerical simulations underestimate the fiber elongation by several orders of magnitude, cf. [33, 30, 25, 36]. While experimental studies show fiber elongations $e\sim\mathcal{O}(10^{6})$, $e=A_{in}/A=d_{in}^{2}/d^{2}$, meaning a reduction of $10^{3}$ in diameter $d$ and of $10^{6}$ in cross-sectional area $A$ compared to the values at the nozzle (indicated by the index in), simulated elongations are of order $e\sim\mathcal{O}(10^{4})$. This is likely due to steady considerations and the neglect of turbulent aerodynamic effects [7, 8, 25]. Assuming an incompressible steady fiber jet the relation $uA=u_{in}A_{in}$ with scalar jet speed $u$ holds true. Hence, the computed elongation is restricted by the velocity $\mathbf{v}_{\star}$ of the surrounding air stream, i.e., $e=u/u_{in}<\lVert\mathbf{v}_{\star}\rVert_{\infty}/u_{in}$. This estimate turns out to be valid also in (instationary) melt-blowing simulations where the surrounding airflow is computed (even) on basis of a turbulence model when only mean airflow informations are taken into account in the aerodynamic driving of the fiber jet [29, 36]. Experiments in [26, 33] indicate the relevance of the turbulent effects for the jet thinning. In [34] a viscoelastic fiber model based on an upper convected Maxwell description (UCM) has been employed for melt-blowing, which is opposed to random pulsations. This is done by applying perturbation frequencies on a rectilinear fiber jet leading to bending instabilities and causing significant stretching and thinning of the jet. The examination has been extended to multiple fibers, focusing on the prediction of fiber deposition patterns and fiber-size distributions in the resulting nonwovens in [35]. Latest works deal with the numerical investigation of the angular fiber distribution, the effect of uptake velocity as well as the lay-down on a rotating drum [27, 14, 13]. In [16] the significance of turbulence for melt-blowing has been approached by studying the effect of turbulent aerodynamic velocity fluctuations on a simplified fiber model of ordinary differential equations. There, a $k$-$\epsilon$ turbulence description of the high-speed airflow serves as basis for the reconstruction of the velocity fluctuations, yielding a stochastic aerodynamic force acting on the fiber jet.

The aim of this paper is to establish a numerical framework for fibers in turbulent air that makes the simulation of industrial melt-blowing processes feasible. For this purpose we bring together the two described approaches: we extend the random field sampling of [16] to the instationary viscoelastic UCM fiber model of [34]. Since the aerodynamic forces are the key player for the fiber behavior, we employ a one-way coupling of the outer air stream with the fibers by the help of the force model given in [22]. Of importance is the efficient and robust realization that enables us presenting numerical results of an industrial setup with an appropriate viscoelastic description of the fiber, the inclusion of temperature effects and the direct incorporation of the turbulence structure of the outer air stream for the first time in literature.

Regarding the viscoelastic UCM fiber model of [34], that is asymptotically derived by slender body theory in [20], in Lagrange description it can uniquely be written as quasilinear hyperbolic first order system of partial differential equations on a growing space-time domain. Its classification with respect to the growing fiber domain gives requirements on boundary conditions with regard to well-posedness of the mathematical problem formulation and suggests a parameterization of the fiber tangent by the help of spherical coordinates. The effects of turbulent fluctuations are calculated by the turbulence reconstruction procedure described in [16] and coupled into the fiber model by an air force function. The resulting instationary problem is solved using finite volumes in space with numerical fluxes of Lax-Friedrichs type as well as employing the implicit Euler method in time. For an industrial melt-blowing setup we show the applicability of our model and numerical solution framework and demonstrate the relevance of the turbulent fluctuations causing fiber elongations of the expected higher order of magnitude compared to stationary simulations. From the repeated random sampling of fibers in the sense of the Monte Carlo method a distribution of the final fiber diameters is obtained that yields fiber diameters of realistic order of magnitude.

The paper is structured as follows. In Sec. 2 we start with the instationary viscoelastic UCM fiber model, regarding its classification and correct closing with boundary conditions. Furthermore, we give a short survey in reconstructing the turbulent fluctuations of an underlying air stream. After that we discuss our numerical solution framework and the handling of the growing fiber domain in Sec. 3. In Sec. 4 we consider an industrial melt-blowing setup, for which we present simulation results covering the turbulent effects due to the high-speed air stream.

2. Viscoelastic fiber melt-blowing model

For the melt-blowing of a fiber in a turbulent air stream we present an asymptotic instationary viscoelastic UCM fiber model in Lagrangian (material) description. We classify the resulting quasilinear system of partial differential equations of first order and discuss the appropriate closing by boundary conditions. The choice of the boundary conditions suggests a description with respect to a fiber tangent associated basis. The fiber tangent $\boldsymbol{\tau}$ with norm $e=\lVert\boldsymbol{\tau}\rVert$ and direction $\mathbf{t}=\boldsymbol{\tau}/\lVert\boldsymbol{\tau}\rVert$ is in particular parameterized by the help of spherical coordinates. Moreover, we present the models for the aerodynamic force and the heat exchange used in the one-way coupling with the surrounding airflow and introduce the stochastic modeling concept by which the effects of the turbulent aerodynamic velocity fluctuations are incorporated in the fiber system.

2.1. Asymptotic fiber jet model

The extrusion of a fiber jet from a nozzle into an air stream can be seen as an inflow problem with a domain enlarging over time. Let $\Omega=\{(\zeta,t)\in\mathbb{R}^{2}\,\lvert\,\zeta\in\mathcal{Q}(t),\,t\in(0,t_{end}]\}$ be the space-time domain with time-dependent growing space $\mathcal{Q}(t)=(-\zeta_{L}(t),0)$, where $\mathrm{d}/\mathrm{d}t\,\zeta_{L}(t)=v_{in}(t)$, $\zeta_{L}(0)=0$, with $v_{in}$ [m/s] being the (scalar) inflow velocity at the nozzle. In the following we assume a constant inflow velocity, i.e., $v_{in}=const$, yielding $\zeta_{L}(t)=v_{in}\,t$. The fiber jet is represented by a time-dependent curve $\mathbf{r}:\Omega\rightarrow\mathbb{R}^{3}$, where the fiber end corresponds to the material parameter $\zeta=0$ and the material points entering the fiber (flow) domain at the nozzle are $\zeta=-\zeta_{L}(t)$. We assume incompressibility of the fiber jet such that besides the mass also the volume is conserved, i.e., $\partial_{t}\varrho_{M}=0$ and $\partial_{t}\varrho_{V}=0$ with mass and volume line densities, $\varrho_{M}$ [kg/m] and $\varrho_{V}$ [m2], respectively. The mass and volume line densities are considered to be constant at the nozzle yielding $\varrho_{M}=\rho d_{in}^{2}\pi/4$ and $\varrho_{V}=\lVert\boldsymbol{\tau}_{in}\rVert d_{in}^{2}\pi/4$ with constant fiber density $\rho$ [kg/m3], nozzle diameter $d_{in}$ [m] and fiber tangent at the nozzle $\boldsymbol{\tau}_{in}$. According to [20], where the viscoelastic UCM string model has been systematically derived by slender body asymptotics, our model for the extruding fiber jet is given in Lagrangian description by

[TABLE]

supplemented with appropriate initial and boundary conditions to be specified. The diameter function $d:\mathcal{Q}\rightarrow\mathbb{R}^{+}$ is introduced via

[TABLE]

The two kinematic equations relate the fiber velocity $\mathbf{v}$ [m/s] and the fiber tangent $\boldsymbol{\tau}$ to the derivatives of the fiber curve $\mathbf{r}$ [m] with respect to time $t$ [s] and material parameter $\zeta$ [m]. The two dynamic equations prescribing the conservation of linear momentum and energy yield equations for fiber velocity $\mathbf{v}$ and fiber temperature $T$ [K]. The acting outer line force densities arise from gravity $\mathbf{f}_{g}=\varrho_{M}g\mathbf{e}_{g}$ [N/m] with direction $\mathbf{e}_{g}$, $\|\mathbf{e}_{g}\|=1$, and gravitational constant $g$ [m/s2] as well as from the surrounding airflow $\mathbf{f}_{air}$ [N/m]. Moreover, $\alpha$ [W/(m2K)] is the heat transfer coefficient, $T_{\star}$ [K] the aerodynamic temperature field, $c_{p}$ [J/(kgK)] the constant specific heat capacity of the fiber and $d$ [m] the fiber diameter. The models for the aerodynamic line force density $\mathbf{f}_{air}$ and for the heat transfer coefficient $\alpha$ are presented in Sec. 2.3. Concerning the viscoelastic material laws, they are based on a UCM model for the fiber stress $\sigma$ [Pa] and pressure $p$ [Pa]. Here, $\mu$ [Pa s] describes the dynamic viscosity and $\theta$ [s] the relaxation time of the fiber jet. Under the assumption of incompressibility the relation $\mu/\theta=E/3$ with elastic modulus $E$ [Pa] holds. We model the dynamic viscosity and relaxation time dependent on the temperature $T$, i.e., $\mu=\mu(T)$, $\theta=\theta(T)$. The corresponding rheological laws for an industrial example are specified in Sec. 4.1.

For the numerical treatment of the problem it is convenient to deal with dimensionless model equations. We introduce the dimensionless quantities as $\tilde{y}(\tilde{\zeta},\tilde{t})=y(\zeta_{0}\tilde{\zeta},t_{0}\tilde{t})/y_{0}$ and use the reference values $y_{0}$ as given in Tab. 2.1. Here, $y_{in}$ indicates the value of a quantity $y$ at the nozzle and $H$ denotes the height of the considered melt-blowing device. The constant mass and volume line densities $\varrho_{M}$, $\varrho_{V}$ become $\tilde{\varrho}_{M}=\tilde{\varrho}_{V}=1$ in dimensionless form. To keep the notation simple we suppress the label $\tilde{~{}}$ in the following. Then, the dimensionless model equations read

[TABLE]

The fiber behavior is characterized by the dimensionless parameters given in Tab. 2.1, that are the Reynolds number $\mathrm{Re}$ as ratio of inertial to viscous forces, the Deborah number $\mathrm{De}$ as ratio of relaxation time to characteristic time-scale, the Froude number $\mathrm{Fr}$ as ratio of inertial to gravitational forces, the Stanton number $\mathrm{St}$ as ratio of heat transfer to thermal capacity as well as the slenderness ratio $\varepsilon$. The time-dependent space domain simplifies to $\mathcal{Q}(t)=(-t,0)$.

Remark 1**.**

The viscoelastic UCM fiber model (2.1) covers the limit cases describing pure viscous as well as elastic material behavior. The limit $\mathrm{De}\rightarrow 0$ yields a viscous fiber model, whereas the limit $\mathrm{Re}\rightarrow 0$, $\mathrm{De}\rightarrow\infty$ with $\mathrm{Re}\mathrm{De}=\mathrm{Ma}^{2}$ describes an elastic behavior. Here, the dimensionless Mach number $\mathrm{Ma}$ is the ratio of inertial to compressive forces.

Remark 2**.**

As pointed out in [20], the pressure $p$ is at least one order of magnitude smaller than the stress $\sigma$ for fibers with high strain rates $\upsilon=\partial_{t}\lVert\boldsymbol{\tau}\rVert/\lVert\boldsymbol{\tau}\rVert\geq 0$ and large Deborah numbers $\mathrm{De}$, in particular $\lvert p\rvert\leq 0.1\sigma$ if $\upsilon\mathrm{De}\,\theta\geq 0.35$. This means the pressure equation can be neglected in such cases. In [34] this simplification is employed instantaneously to the UCM model for melt-blowing.

2.2. Classification and boundary conditions

The dimensionless fiber model (2.1) can uniquely be written as a quasilinear system of partial differential equations of first order [20]

[TABLE]

with the vector of unknowns $\boldsymbol{\varphi}=(\mathbf{r},\boldsymbol{\tau},\mathbf{v},T,\sigma,p)\in\mathbb{R}^{12}$. The system is classified mathematically by the spectrum of the system matrix $\mathbf{M}$ that consists of the eigenvalues

• •

$\lambda_{1}=0$ (multiplicity $6$),

• •

$\lambda_{2,3}=\pm\sqrt{\sigma}/\lVert\boldsymbol{\tau}\rVert$ (multiplicity $2$ each),

• •

$\lambda_{4,5}=\pm\sqrt{w}/\lVert\boldsymbol{\tau}\rVert$ (multiplicity $1$ each), $\qquad w=\left(3\mu/\theta+\mathrm{Ma}^{2}\left(\sigma+3p\right)\right)/\mathrm{Ma}^{2}$

The system is of hyperbolic type if $\sigma>0$ and $w>0$. Otherwise, it is mixed elliptic-hyperbolic, or even shows a parabolic deficiency if $\sigma=0$ and/or $w=0$.

Since the hyperbolic case is relevant for the application, we focus on it and discuss the closing of the system by appropriate boundary and initial conditions. At the fiber jet end, which corresponds to a fixed material point in Lagrangian description ($\zeta=0$), the characteristic related to the eigenvalue $\lambda_{i}$ runs from the nozzle to the jet end if $\lambda_{i}>0$ and from the jet end towards the nozzle if $\lambda_{i}<0$. At the nozzle ($\zeta=-\zeta_{L}(t)$) the orientations of the characteristics depend on the scalar inflow velocity of the fiber jet, which reads $v_{in}/v_{0}=1$ in non-dimensional form. If $\lambda_{i}>-v_{in}/v_{0}=-1$ for $i\in\{1,...,5\}$, the corresponding characteristic propagates from the nozzle to the jet end, otherwise the other way round. The orientations of the characteristics yield requirements on the boundary conditions with regard to the well-posedness of the problem. Since $\lambda_{3}<0$ (multiplicity 2) and $\lambda_{5}<0$ (multiplicity 1), we have to pose three boundary conditions at the fiber jet end. Because of the spinning setup we model the fiber end ($\zeta=0$) as stress-free, i.e.,

[TABLE]

Employing the viscoelastic material law for $\sigma$ yields a constant fiber elongation $e=\lVert\boldsymbol{\tau}\rVert$ at the fiber end over time, i.e., $\partial_{t}e(0,t)=0$. To preserve this compatibility condition we pose

[TABLE]

assuming the fiber jet to leave the nozzle unstretched. The eigenvalues $\lambda_{1}$, $\lambda_{2}$, $\lambda_{4}$ are non-negative and thus imply nine boundary conditions at the nozzle ($\zeta=-\zeta_{L}(t)$, $t\geq 0$),

[TABLE]

Here, $\mathbf{r}_{in}$ is assumed to be constant. Furthermore, we set the following initial conditions for $t=0$,

[TABLE]

Depending on the propagation-speed of the characteristics at the nozzle we pose further boundary conditions: we additionally prescribe for $t>0$

[TABLE]

The total time-derivative of the fiber curve $\mathbf{r}$ at the nozzle yields the compatibility condition $\mathbf{v}(-\zeta_{L}(t),t)=\boldsymbol{\tau}(-\zeta_{L}(t),t)$ for all times $t$. Through the above choice of the boundary conditions for $\mathbf{v}$, $\boldsymbol{\tau}/e$, and $e$ at the nozzle this condition is inherently fulfilled.

The choice of the boundary conditions and in particular the decomposition of the fiber tangent $\boldsymbol{\tau}=e\mathbf{t}$ into elongation $e=\lVert\boldsymbol{\tau}\rVert$ and direction $\mathbf{t}$, $\|\mathbf{t}\|=1$, suggests a reformulation of the corresponding dynamic equation $\partial_{\zeta}\mathbf{r}=\boldsymbol{\tau}$. Making use of the compatibility condition $\partial_{t}\boldsymbol{\tau}=\partial_{t}\partial_{\zeta}\mathbf{r}=\partial_{\zeta}\partial_{t}\mathbf{r}=\partial_{\zeta}\mathbf{v}$ yields an equation for the elongation $e$

[TABLE]

The normalized tangent $\mathbf{t}$ can be parameterized by means of spherical coordinates

[TABLE]

Then, its time-derivative reads $\partial_{t}\mathbf{t}=\mathbf{n}\partial_{t}\vartheta+\mathbf{b}\partial_{t}\varphi$ with normal $\mathbf{n}=(\cos\vartheta\cos\varphi,\cos\vartheta\sin\varphi,-\sin\vartheta)$ and binormal $\mathbf{b}=(-\sin\vartheta\sin\varphi,\sin\vartheta\cos\varphi,0)$. The set $\{\mathbf{t},\mathbf{n},\mathbf{b}\}\subset\mathbb{R}^{3}$ forms an orthogonal basis where $\|\mathbf{t}\|=\|\mathbf{n}\|=1$. Employing

[TABLE]

gives relations for the polar $\vartheta$ and azimuth angles $\varphi$

[TABLE]

Summing up, our viscoelastic instationary fiber model on a growing domain in Lagrangian description is given by System 3.

System 3** (Instationary viscoelastic fiber model).**

Kinematic and dynamic equations as well as material laws in $\Omega$:

[TABLE]

Initial-boundary conditions at the nozzle ($\zeta=-\zeta_{L}(t)$, $t\geq 0$):

[TABLE]

Initial conditions ($t=0$):

[TABLE]

Boundary conditions at the nozzle ($\zeta=-\zeta_{L}(t)$, $t>0$):

[TABLE]

Boundary conditions at the fiber end ($\zeta=0$, $t>0$):

[TABLE]

2.3. Exchange models for one-way coupling with turbulent airflow

In this work we consider a one-sided coupling of the airflow with the fiber, neglecting feedback effects of the fiber on the airflow. The respective exchange models used for the aerodynamic line force density $\mathbf{f}_{air}$ and the heat transfer coefficient $\alpha$ are briefly summarized in this subsection. Moreover, we describe the concept how the turbulent aerodynamic velocity fluctuations are realized with respect to an underlying (stochastic) airflow simulation and incorporated in our fiber model (System 3).

Note that to distinguish the fiber quantities from the airflow quantities, all airflow associated fields are labeled with the index ⋆ as before. In particular, $\mathbf{v}_{\star}$ denotes the velocity, $\rho_{\star}$ the density, $\nu_{\star}$ the kinematic viscosity, $c_{p,\star}$ the specific heat capacity, $\lambda_{\star}$ the thermal conductivity, $k_{\star}$ the turbulent kinetic energy and $\epsilon_{\star}$ the viscous dissipation of the turbulent motions per unit mass of the air. 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 Tab. 2.1.

2.3.1. Aerodynamic force and heat transfer coefficient

The models for the aerodynamic force and the heat transfer coefficient are determined by material and geometrical properties as well as the incident flow situation which can be prescribed by the fiber orientation (normalized tangent) $\mathbf{t}$ and the relative velocity between airflow and fiber $\mathbf{v}_{\star}-\mathbf{v}$.

The aerodynamic line force density $\mathbf{f}_{air}$ is modeled by means of a dimensionless drag function $\mathbf{F}:\mathrm{SO}(3)\times\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$ which depends on fiber tangent and relative velocity,

[TABLE]

The drag function can be particularly expressed in terms of the tangential $\mathbf{w_{t}}=(\mathbf{w}\cdot\mathbf{t})\mathbf{t}$ and normal relative velocity components $\mathbf{w_{n}}=\mathbf{w}-\mathbf{w_{t}}$, $w_{n}=\lVert\mathbf{w_{n}}\rVert$. The models used for the tangential and normal air resistance coefficients $r_{t}$, $r_{n}$ are taken from [22], see Appendix A.1 for details. The occurring dimensionless numbers are the air drag associated number $\mathrm{A}_{\star}$ and the mixed (air-fiber) Reynolds number $\mathrm{Re}_{\star}$ (cf. Tab. 2.1). Concerning lift forces see Remark 4.

The heat transfer coefficient $\alpha$ is modeled by a Nusselt number associated dimensionless function $\mathcal{N}:\mathbb{R}^{3}\rightarrow\mathbb{R}$ which depends on the tangential and absolute relative velocity and the Prandtl number,

[TABLE]

For details on the used heuristic model for $\mathcal{N}$ we refer to Appendix A.2. The occurring dimensionless numbers are the Nusselt number $\mathrm{Nu}_{\star}$, the Prandtl number $\mathrm{Pr}_{\star}$ as well as the mixed (air-fiber) Reynolds number $\mathrm{Re}_{\star}$ (cf. Tab. 2.1).

2.3.2. Turbulence reconstruction

A direct numerical simulation of the turbulent airflow in the application is not possible due to the required high resolution. Hence, a statistical turbulence description is used where the airflow velocity $\mathbf{v}_{\star}$ is assumed to consist of a mean (deterministic) part $\bar{\mathbf{v}}_{\star}$ and a fluctuating (stochastic) part $\mathbf{v}^{\prime}_{\star}$, i.e.,

[TABLE]

The mean velocity is given by the Reynolds-averaged Navier-Stokes equations, while the fluctuations are only characterized by certain quantities that the respective turbulence model provides. To obtain $\mathbf{v}^{\prime}_{\star}$ explicitly as random field we apply a turbulence reconstruction that has been developed in [16] on basis of a $k_{\star}$-$\epsilon_{\star}$ turbulence model. Assuming given dimensionless space-time-dependent fields for the turbulent kinetic energy $k_{\star}$ and the viscous dissipation of the turbulent motions per unit mass $\epsilon_{\star}$, the general concept of the turbulence reconstruction is to model the local turbulent fluctuations as homogeneous, isotropic, incompressible Gaussian random fields in space and time, $\mathbf{v}_{\star,loc}^{\prime}=\mathbf{v}_{\star,loc}^{\prime}(\mathbf{x},t;\nu_{\star},\bar{\mathbf{v}}_{\star})$, that depend parametrically on the kinematic viscosity and mean velocity of the airflow, as done in [21, 22]. To form the large-scale structure of the global turbulence the local fluctuations fields are superposed based on a Global-from-Local assumption. The globalization strategy according to [16] yields

[TABLE]

Besides the slenderness ratio $\varepsilon$ and the mixed Reynolds number $\mathrm{Re}_{\star}$, the occurring dimensionless numbers are the degree of turbulence $\mathrm{Tu}_{\star}$ and the turbulent time scale ratio $\mathrm{Tt}_{\star}$ as given in Tab. 2.1.

Note that the occurring turbulent length and time scales give requirements on the spatial and temporal resolution in our numerical solution algorithm (cf. Rem. 7 in Sec. 3). In particular, $l^{\prime}_{\star}=\mathrm{Tu}_{\star}/\mathrm{Tt}_{\star}k_{\star}^{3/2}/\epsilon_{\star}$ is the dimensionless turbulent length scale indicating the expected length of the large-scale vortices, and $t_{\star}^{\prime}=1/\mathrm{Tt}_{\star}k_{\star}/\epsilon_{\star}$ is the dimensionless turbulent time scale describing the expected creation and break-up time of the vortices. For details on the general sampling procedure providing a fast and accurate sampling of the random fields we refer to [16]. To even increase the efficiency of the procedure we use here a simplified underlying energy spectrum, see Appendix B for details on the modeling of $\mathbf{v}_{\star,loc}^{\prime}$.

Remark 4** (Lift forces).**

In industrial melt-blowing processes lift forces on a fiber are created through airflow vortices approaching the fiber and by vortex shedding at the back of the fiber. While the latter can be neglected since the fiber is meanly following the turbulent air stream, the first mechanism is included by the help of the following ansatz: the local turbulent instationary velocity fluctuations $\mathbf{v}_{\star}^{\prime}$ are plugged into the air drag model (2.3), meaning local observations are mapped into a stationary far field consideration. This leads to aerodynamic forces on the fiber acting perpendicular to the $(\bar{\mathbf{v}}_{\star}-\mathbf{v})$-$\boldsymbol{\tau}$-plane.

3. Numerical Scheme

System 3 is a boundary value problem of a quasilinear system of partial differential equations of first order on a growing domain. It is discretized with finite volumes in space based on a central flux approximation with a Lax-Friedrich type stabilization and with the implicit Euler method in time. The growing fiber domain is realized by dynamic and static spatial cells according to the discretization concept in [1].

We reformulate System 3 as

[TABLE]

with the vector of unknowns $\mathbf{y}=(\mathbf{r},e,\vartheta,\varphi,\mathbf{v},T,\sigma,p)\in\mathbb{R}^{12}$ and consider it on the spatial domain $\mathcal{Q}(t)=(-t,0)$ for times $0\leq t\leq t_{end}$. The introduction of the matrix $\mathbf{K}$ avoids a singularity for $\sin\vartheta=0$. For $\sin\vartheta\neq 0$, $\mathbf{K}$ is invertible revealing the unique quasilinear form (2.2).

For the spatial discretization we employ a finite volume scheme. We introduce a constant cell size $\Delta\zeta$ and define the number of dynamic cells $N(t)$ depending on the fiber length $\zeta_{L}(t)=t$ at time $t$ as

[TABLE]

where $\lfloor\cdot\rfloor$ denotes the floor function. Furthermore, we introduce the discretization points

[TABLE]

The points $\zeta_{i}$, $i=1,...,N(t)$, represent the cell centers. The dynamic cell closest to the nozzle ($\zeta=-t$) is given by $[\zeta_{1/2},\zeta_{3/2}]$, whereas $\zeta_{N+1/2}=0$ is the fiber end, cf. Fig. 3.1. The jet growth is realized by adding static cells at the nozzle. This means we add new cells, which are initialized by the boundary conditions at the nozzle (i.e., at the left side of the computational domain). The cells remain static until they have completely left the nozzle. When they have completely entered the flow domain they are called dynamic cells and are then taken into consideration for the computation. The introduction of static cells at the nozzle allows the suitable initialization of a jet with length $\zeta_{L}(t)<\Delta\zeta$ and a stable numerical treatment of the temporal evolution.

We define the cell averages $\mathbf{y}_{i}$, $i=1,...,N(t)$, of the unknown quantities as

[TABLE]

integrate the quasilinear system (3.1) over the control cells $[\zeta_{i-1/2},\zeta_{i+1/2}]$, $i=1,...,N(t)$, in which we assume $\mathbf{X}(\mathbf{y})|_{[\zeta_{i-1/2},\zeta_{i+1/2}]}=\mathbf{X}(\mathbf{y}_{i})$ for $\mathbf{X}=\mathbf{K},\mathbf{L},\mathbf{l}$ and adopt the idea of the Lax-Friedrichs scheme for the approximation of the numerical fluxes as done in [12]. The resulting system of ordinary differential equations for the cell averages $\mathbf{y}_{i}$ with respect to time has the form

[TABLE]

where $\Delta t$ denotes the constant time-step size, that we will use in the temporal discretization. The incorporation of initial-boundary and boundary conditions in our numerical scheme is realized by ghostlayers. Following [18] quantities not being prescribed at a boundary are extrapolated on the corresponding ghostlayer, in particular we choose first order extrapolation.

For the solution of the system of ordinary differential equations (3.2) we employ the stiffly accurate implicit Euler scheme with constant time-step size $\Delta t$

[TABLE]

with $\mathbf{y}_{i}^{n}=\mathbf{y}_{i}(t^{n})$ and $t^{n}=n\Delta t$ for $n=0,...,M$, $t^{M}=t_{end}$. The resulting nonlinear system of equations is solved by a Newton-method with Armijo step size control, where the Jacobian of the system matrix is prescribed analytically. The break-up criterion of the Newton algorithm is set to an absolute error tolerance $tol=10^{-8}$ with respect to the maximum norm.

Remark 5** (Artificial diffusion).**

The semi-discrete system (3.2) can be seen as a spatial discretization of

[TABLE]

with $\eta=(\Delta\zeta)^{2}/(2\Delta t)$ by means of a central approximation of the flux terms. This means we add artificial diffusion of magnitude $\eta$ to our system as it is common for classical Lax-Friedrich schemes.

Remark 6** (Convergence of numerical scheme).**

As it is well-known from hyperbolic literature (e.g. [18]), the numerical scheme (3.3) provides accuracy of order one with respect to the time and accuracy of order two with respect to the space discretization yielding a combined convergence of order one. In [9, 12, 23] a similar scheme has been investigated with respect to a stability concept for non-conservative hyperbolic partial differential equations.

Remark 7** (Spatial and temporal resolution).**

The temporal and spatial grid sizes have to be chosen in such a way that the turbulent scales of the underlying airflow are resolved properly. In particular, the turbulent length scale $l^{\prime}_{\star}=\mathrm{Tu}_{\star}/\mathrm{Tt}_{\star}k_{\star}^{3/2}/\epsilon_{\star}$ and the turbulent time scale $t_{\star}^{\prime}=1/\mathrm{Tt}_{\star}k_{\star}/\epsilon_{\star}$ used in the turbulence reconstruction (2.5) have to be considered. Furthermore, the time that a vortex needs to pass a fixed material point of the fiber due to their relative velocity has to be taken into account for the temporal resolution. In total, the requirements for a successful simulation in terms of $\Delta\zeta$ and $\Delta t$ read

[TABLE]

Appropriate grid sizes are estimated by computing the bounds for all given airflow data with assumptions on the maximal fiber velocity and elongation.

4. Industrial melt-blowing simulation

In this section we investigate an industrial melt-blowing scenario that has been studied in [16] by the help of a simplified ODE model for the fiber jet position, velocity and elongation. We employ our more sophisticated PDE fiber jet model (System 3), which additionally contains a viscoelastic material behavior and thermal effects describing the jet cooling and solidification. Before we present our simulation results, we specify the industrial setup and state the closing models for the dynamic viscosity as well as relaxation time and elastic modulus. In the scenario we face step size restrictions (cf. Remark 7) that prevent the computability of the whole fiber from nozzle to conveyor belt. To handle this numerical problem we suggest and discuss an appropriate simulation strategy.

4.1. Setup and model closing

In the melt-blowing setup a high-speed air stream is directed vertically downwards in direction of gravity and enters the domain of interest via thin slot dies. The spinning nozzles are located in between and extrude the polymeric fiber jets in the same direction, see Fig. 4.1. We choose an outer orthonormal basis $\{\mathbf{e}_{x},\mathbf{e}_{y},\mathbf{e}_{z}\}$, where $\mathbf{e}_{x}$ points against the direction of gravity (i.e., $\mathbf{e}_{x}=-\mathbf{e}_{g}$) and $\mathbf{e}_{y}$ is aligned with the slot inlet. The mean quantities of the turbulent airflow are time-independent and homogeneous in $\mathbf{e}_{y}$-direction, such that a stationary $k_{\star}$-$\epsilon_{\star}$ simulation for a representative two-dimensional cut showing the $\mathbf{e}_{x}$-$\mathbf{e}_{z}$-plane is reasonable. The origin of the outer basis is aligned with the external given airflow data, such that the considered nozzle is at the position $\mathbf{r}_{in}=(-2.85\cdot 10^{-2},0,0)\,$m in the airflow field. We use the same $k_{\star}$-$\epsilon_{\star}$ simulation results as in [16], supplemented with an additional temperature profile as depicted in Fig. 4.2. The melt-blown fiber polymer is of polypropylene (PP) type with material parameters taken from [16]. The process and physical parameters as well as the resulting dimensionless numbers are listed in Table 4.1.

For the temperature-dependent dynamic viscosity of the PP-type fiber material we employ the Arrhenius law. The corresponding relation for the dimensionless viscosity $\mu$ depending on the temperature $T$ is given by

[TABLE]

The polymer-specific constants coming from measurements are $a_{\mu}=0.1352$ Pas, $b_{\mu}=852.323$ K, $c_{\mu}=273.15$ K. We choose the following heuristic model for the relaxation time

[TABLE]

with $a_{\theta}=10^{9}$ Pa and $b_{\theta}=2\cdot 10^{8}$ Pas showing a meaningful limit behavior: for $T\rightarrow\infty$ the dimensional relaxation time is of order $\mathcal{T}\sim\mathcal{O}(10^{-1}$ s), which is typical for melt-blown polymers, see, e.g., [34]. Furthermore, employing the relation $\mu/\theta=E/3$ the resulting dimensionless elastic modulus $E$ reads

[TABLE]

For $T\rightarrow c_{\mu}=273.15$ K the dimensional elastic modulus $\mathcal{E}$ approaches $\mathcal{E}=a_{\theta}=10^{9}$ Pa – a typical value for hardened polymer, see, e.g., [5].

4.2. Simulation strategy

Expecting a maximal fiber elongation $e=10^{6}$ and a maximal dimensionless relative velocity between fiber and airflow $\lVert\mathbf{v}_{\star}-\mathbf{v}\rVert\leq\lVert\mathbf{v}_{\star}\rVert=4.78\cdot 10^{4}$ (in dimensionless form) in the industrial melt-blowing, the step size restriction for the spatial and temporal fiber discretization (3.4) gives

[TABLE]

(cf. Fig. 4.2). Such a resolution implies computationally impractical runtimes. However, to make a simulation for the setup feasible, we suggest the following strategy that is motivated from observations of the process.

In the region close to the nozzle the high-speed air stream pulls the slowly extruded fiber jet rapidly down without any lateral bending. The hot temperatures prevent fiber cool-down and solidification. Thus, the magnitude of the Deborah number $\mathrm{De}$ at the nozzle (cf. Tab. 4.1) allows the consideration of the viscous limit case $\mathrm{De}\rightarrow 0$ (see Remark 1). Moreover, the fiber jet behavior is mainly determined by the mean airflow, turbulent effects play a negligible role. Hence, we assume that in the nozzle region (i.e., deterministic region) the polymer jet can be described by a steady uni-axial viscous fiber model with deterministic aerodynamic force and heat transfer. This model follows from System 3 by a re-parameterization into Euler (spatial) description, transition to steady-state, and the limit $\mathrm{De}\rightarrow 0$. The resulting boundary value problem of ordinary differential equations is solved by a continuation-collocation method, which we successfully employed in studies on glass wool manufacturing [3], electrospinning [2] and dry spinning [31]. Further details on the model and its numerical treatment are given in Appendix C. Note that the use of the viscous fiber model is not only physically reasonable, but it also simplifies crucially the numerical treatment. Concerning the viscoelastic fiber model, the rapid changes of the fiber quantities in the nozzle region caused by the immediate pull down of the fiber yields multiple changes in the structure of the quasilinear system matrix in view of its eigenvalues and its resulting classification. This means that the runs of the characteristics change their direction several times. In the steady uni-axial model this leads to singular system matrices and closing problems with appropriate boundary conditions making the numerical treatment extremely complicated. This issue has been addressed by [19] in the context of existence regimes for solutions of an uni-axial UCM fiber model under gravitational forces. We circumvent these problems when using the viscous fiber model where no mathematical regime changes take place.

In the region away from the nozzle the turbulent aerodynamic fluctuations crucially affect the fiber behavior (i.e., stochastic region). By means of the uni-axial steady fiber solution (from the nozzle region) we identify a coupling point, from where on the further fiber behavior downwards to the bottom is described by the instationary viscoelastic fiber model (System 3) accounting for turbulent effects. The simulation with the numerical scheme from Sec. 3 becomes here feasible since the expected elongation and relative velocities in this domain are much smaller and hence the spatial and temporal step size restrictions weaken compared to (4.1).

The coupling between the stationary and instationary fiber simulations is done in the following way: Let the fiber domain in the Eulerian parametrization $\Omega(t)=\Omega_{d}\cup\Omega_{s}(t)$ be divided into the time-independent deterministic part $\Omega_{d}$, where the fiber is uni-axially stretched, and the time-dependent stochastic part $\Omega_{s}(t)$, where the fiber is strongly affected by the turbulent fluctuations. Consider $C=\Omega_{d}\cap\Omega_{s}(t)$ to be the time-independent coupling point between the deterministic and the stochastic domain. First, we perform the simulation of the steady viscous fiber model (System 8 in Appendix C) for the whole fiber domain, i.e., $\Omega_{d}=\Omega$, yielding solutions for the scalar fiber speed $u$, temperature $T$, stress $\sigma$ and pressure $p$. Second, we determine the coupling point $C$ by the ratio of the relative velocity between fiber and airflow $v_{rel}=\lVert\mathbf{v}_{\star}-u\boldsymbol{\tau}\rVert$ and the turbulent velocity scale $k_{\star}^{1/2}$, in particular

[TABLE]

So the coupling point is the nearest point to the nozzle, where the ratio of the relative velocity and the turbulent velocity scale is below one order of magnitude. At $C$ the quantities of the stationary solution are denoted by $u_{C}$, $T_{C}$, $\sigma_{C}$, $p_{C}$. Third, for the subsequent solving of the instationary viscoelastic fiber model (System 3) on $\Omega_{s}(t)$ (reformulated in Lagrangian coordinates) we adjust the typical values and adapt the initial conditions. We particularly set the reference values used for the non-dimensionalization (see Sec. 2.1) to be

[TABLE]

then the dimensionless numbers change accordingly. The altered initial conditions read

[TABLE]

These modifications can be interpreted as putting a fictive nozzle with adjusted extrusion conditions at the spatial position of the coupling point $C$. The diameter of the fictive nozzle reflects the pre-elongations of the extruded fiber by the factor $v_{C}/v_{in}$ compared to $d_{in}$ of the original nozzle.

In the setup the crucial stretching of the fiber takes place in the upper part of the device and ends when the fiber is nearly solidificated. Since we are interested in the maximal achieved fiber elongations as well as in the corresponding fiber diameter distribution, but not in the lay-down process, it is sufficient to cutoff the fiber before it reaches the bottom of the device. We choose to cutoff the fiber, when it reaches the height corresponding to $x=-9.45\cdot 10^{-2}$ m. Below this point the airflow temperature satisfies $T_{\star}<353.15$ K (see Fig. 4.2). We expect the fiber dynamic viscosities to be of magnitude $\mu\sim\mathcal{O}(10^{3}$ Pas$)$, implying that no noticeable further fiber elongations take place.

4.3. Results

In the following we present the numerical results for the industrial spinning setup described in Sec. 4.1. While all computations for a single fiber realization have been done on an Intel Core i7-6700 CPU (4 cores, 8 threads) and 16 GBytes of RAM, the Monte Carlo simulation has been performed on a MPI cluster (dual Intel Xeon E5-2670, 16 CPU cores per node, 64 GB RAM) with one CPU core for each fiber realization. For all computations the MATLAB version R2016b has been used.

Figure 4.3 shows the results for scalar speed $u$, stress $\sigma$, pressure $p$ and temperature $T$ as well as the induced elongation $e$ of the steady viscous uni-axial fiber model in an Eulerian parameterization that spans the whole domain $\Omega$. The maximal fiber speed is $u=1.74\cdot 10^{2}$ m/s and the corresponding maximal fiber elongation is $e=1.74\cdot 10^{4}$. This indicates that a stationary fiber simulation is not physically reasonable for the whole domain $\Omega$, since much higher fiber elongations for the melt-blowing setup are expected. Nevertheless, the steady viscous solution serves as adequate approximation of the fiber behavior in the nozzle region as described in Sec. 4.2. For the further instationary viscoelastic simulation we determine the spatial position of the coupling point $C$ and put a fictive nozzle at $\mathbf{r}_{in}=(-3.83\cdot 10^{-2},0,0)\,$m. The corresponding fiber quantities at this fictive nozzle are

[TABLE]

and the dimensionless numbers change accordingly, see Tab. 4.2.

Considering the stochastic region, the numerical step size restriction (3.4) for the fiber discretization weakens compared to (4.1)

[TABLE]

we choose $\Delta\zeta=\Delta t=10^{-5}$ for our computation. As expected the turbulent fluctuations of the airflow cause a swirling of the fiber jet such that the fiber curve leaves the $\mathbf{e}_{x}$-axis shortly away from the fictive nozzle. Figure 4.4 shows temporal snapshots of the curve for one representative fiber before its cutoff (at $x=-9.45\cdot 10^{-2}$ m). The fluctuations move the fiber jet not only downwards but also upwards such that the fiber curve creates loops. In these loops high aerodynamic forces act on the fiber due to high relative velocity gradients causing the fiber to elongate. Figure 4.5 shows exemplary the temporal evolution of the fiber quantities for one material point. Obviously the material point experiences high elongations: directly after entering the flow domain high relative velocities in $\mathbf{e}_{x}$-direction between the fiber velocity $v_{1}$ and the deterministic airflow velocity $\bar{v}_{\star,1}$ cause a fiber stretching. After the fiber velocity $v_{1}$ reaches the corresponding deterministic airflow velocity $\bar{v}_{\star,1}$ the fiber experiences a further stretching due to the velocity fluctuations, in which the mean stretching takes place in regions where high lateral air velocities $v_{\star,2}$, $v_{\star,3}$ create swirls. The final elongation at this material point is of magnitude $e\sim\mathcal{O}(10^{5})$ and therewith clearly exceeds the theoretically possible deterministic expectations. In particular, the computed elongation in a stationary simulation is obviously restricted by the velocity of the air stream, i.e., $e_{max}=u/v_{in}<\lVert\mathbf{v}_{\star}\rVert_{\infty}/v_{in}=4.78\cdot 10^{4}$. Furthermore, the stationary uni-axial viscous simulation only achieves $e=1.74\cdot 10^{4}$ (cf. Fig. 4.3). In the region of high fiber stretching the material point experiences high stresses $\sigma$ that partly dissipate due to the elastic material behavior before the fiber completely solidificates. The pressure $p$ is orders of magnitude smaller compared to the stress $\sigma$ and could therefore be neglected in the simulation as already pointed out in Remark 2. The fiber temperature $T$ approaches the air temperature $T_{\star}$ leading to a cool-down and induced solidification of the jet.

When the fiber reach the cutoff height $x=-9.45\cdot 10^{-2}$ m at time $t=3.92\cdot 10^{-4}$ s, we cutoff the fiber end, track the fiber elongations $e$ as well as the corresponding fiber diameters $d$ and document the occurring relative frequencies until the end time $t_{end}=2.20\cdot 10^{-2}$ s is reached, see Fig. 4.6. To achieve comparability with experiments, we weight the relative frequencies of the fiber diameters with the associated fiber elongations $e$ leading to a diameter distribution in the sense of an Eulerian (spatial) parameterization of the fibers. The resulting elongation and fiber diameter distributions are computed by the help of a Monte Carlo simulation with $67$ samples. We observe a mean elongation $e=9.47\cdot 10^{4}$ again exceeding the deterministic expectations. The mean fiber diameter is $d=1.28\cdot 10^{-6}$ m. This is a typical value for fibers produced in industrial melt-blowing setups, see for example [11]. So our instationary viscoelastic fiber model using an adjusted nozzle as well as employing fluctuation reconstruction of the underlying turbulence effects from an airflow simulation predicts quantitatively well the fiber jet thinning observed in experiments, which would not be possible with a pure steady deterministic simulation neglecting the turbulent aerodynamics velocity fluctuations.

Summing up, our proposed procedure makes the simulation of industrial melt-blowing processes with inclusion of turbulent and viscoelastic effects as well as temperature dependencies feasible. Including turbulent effects acting on the fiber by the help of reconstructing the turbulent structure of the outer air stream yields a jet thinning exceeding the deterministic expectations and produces final fiber diameters of realistic order of magnitude. So our presented modeling and solution framework provides the basis for further parameter studies and the optimization of melt-blowing processes. The computation time for the presented setup is around $96.4$ hours. A combined experimental and numerical study is left to future research.

5. Conclusion

In this paper we presented a model and simulation framework that allowed the numerical investigation of the physical mechanism being responsible for the strong fiber thinning in industrial melt-blowing processes. Considering an asymptotic instationary viscoelastic UCM fiber jet model driven by turbulent aerodynamic forces, the random field sampling strategy of [16] provides an efficient numerical procedure for the realization of the turbulent air flow fluctuctuations. The computed fiber diameters are much lower than those obtained from previous stationary simulations regarding a pure deterministic aerodynamic force on the fiber. Our simulation results clearly stress the significance of the turbulent effects as key player for the production of fibers of micro- and nanoscale. Further parameter studies and an optimization of the industrial process setup will provide the opportunity of simulating fibers with elongations of order $e\sim\mathcal{O}(10^{6})$ compared to the nozzle diameter. In view of more quantitative predictions of the resulting nonwovens a combined experimental and numerical study with experimentally measured temperature-dependencies of polymer properties (e.g., relaxation time) is aimed at in future.

Acknowledgments

This work has been supported by German DFG, project 251706852, MA 4526/2-1, WE 2003/4-1.

Appendix A Exchange models between airflow and fiber

The exchange models between airflow and fiber (Sec. 2.3.1) depend on the (in-)flow situation prescribed by the fiber orientation (normalized tangent) $\mathbf{t}$ and the relative velocity between airflow and fiber $\mathbf{w}$. They go back to studies of a stationary perpendicular laminar flow around a cylinder and have been extended to cover arbitrary angle of attacks and velocity regimes.

A.1. Air resistance coefficients

The dimensionless air drag function $\mathbf{F}:\mathrm{SO}(3)\times\mathbb{R}^{3}\rightarrow\mathbb{R}^{3}$, $\mathbf{F}(\mathbf{t},\mathbf{w})=r_{n}(w_{n})\mathbf{w_{n}}+r_{t}(w_{n})\mathbf{w_{t}}$ can be expressed in terms of the normal $\mathbf{w_{n}}=\mathbf{w}-\mathbf{w_{t}}$, $w_{n}=\lVert\mathbf{w_{n}}\rVert$ and tangential $\mathbf{w_{t}}=(\mathbf{w}\cdot\mathbf{t})\mathbf{t}$ relative velocity components. For the air resistance coefficients $r_{n}$, $r_{t}$ we use the following model taken from [22]

[TABLE]

with $S(w_{n})=2.0022-\log(w_{n})$. It is composed of asymptotic Oseen theory, Taylor heuristic and simulations where the transition points $w_{1}=0.1$ and $w_{2}=100$ are estimated from a least-square approximation of experimental and numerical data. For tangential incident flow situations ($w_{n}\rightarrow 0$) a regularization employs the Stokes theory, yielding the Stokes limits

[TABLE]

as well as the transition point w_{0}=2\exp\big{(}2.0022-{4\pi}/{r_{n}^{S}}\big{)} with the regularization parameter $\epsilon=3.5\cdot 10^{-2}$. The parameters $p_{k,j}$ and $q_{k,j}$ ($k\in\{n,t\}$, $j\in\{0,1,2,3\}$) ensure smoothness,

[TABLE]

[TABLE]

A.2. Heat transfer coefficient

For the heat transfer coefficient $\alpha$ we model the Nusselt number associated function $\mathcal{N}:\mathbb{R}^{3}\rightarrow\mathbb{R}$ as

[TABLE]

depending on the tangential and absolute relative velocity ($w_{t}=\|\mathbf{w_{t}}\|$, $w=\|\mathbf{w}\|$) and the Prandtl number, with $\delta=7.3\cdot 10^{-5}$. The model goes back to [28] where originally a stationary perpendicular laminar flow situation ($w_{t}=0$) around a cylinder for $wp\geq\delta$ was studied and has been extended to ensure a smooth transition to the limit value $\mathcal{N}\rightarrow\mathcal{N}_{0}=0.1$ for vanishing $wp\rightarrow 0$,

[TABLE]

with the coefficients

[TABLE]

and the constant parameters

[TABLE]

The incorporation of the function ${h}$ accounts for varying incident flow directions. It is mainly the squared cosine of the angle of attack, which is regularized to ensure smoothness for tangential incident flow situations with regularization parameter $\epsilon=10^{-7}$,

[TABLE]

Appendix B Model for local turbulent velocity fluctuations

The turbulence reconstruction (cf. Sec. 2.3.2) goes back to the works [21, 22]. Given a $k_{\star}$-$\epsilon_{\star}$ description of the turbulent airflow, the local turbulent velocity fluctuations are modeled as homogeneous isotropic incompressible Gaussian random field in space and time $\mathbf{v}_{\star,loc}^{\prime}$ (with expectation $\mathbb{E}(\mathbf{v}_{\star,loc}^{\prime})=\mathbf{0}$) that depend parametrically on the local kinematic viscosity and mean velocity of the airflow. Its covariance function is prescribed as product of the initial correlations transported with the local mean velocity and their temporal decay

[TABLE]

with $\boldsymbol{\gamma}:\mathbb{R}^{3}\times\mathbb{R}\rightarrow\mathbb{R}^{3}$ isotropic and $\Phi:\mathbb{R}\rightarrow\mathbb{R}$ even for any $\mathbf{x},\mathbf{y},\mathbf{w}\in\mathbb{R}^{3}$ and $t,s,\xi\in\mathbb{R}^{+}_{0}$. Regarding the local kinetic turbulent energy and dissipation rate the following relations (in dimensionless formulation) hold

[TABLE]

where $:$ denotes the scalar product for matrices and $\nabla_{\mathbf{x}}$ the nabla operator with respect to the variable $\mathbf{x}$. Due to the assumptions of isotropy and incompressibility the Fourier-transform $\mathcal{F}_{\boldsymbol{\gamma}}$ of the initial correlations $\boldsymbol{\gamma}$ can be expressed in terms of a single scalar-valued function, the so-called energy spectrum $E:\mathbb{R}^{+}_{0}\times\mathbb{R}^{+}_{0}\rightarrow\mathbb{R}^{+}_{0}$, i.e.,

[TABLE]

with unit matrix $\mathbf{I}\in\mathbb{R}^{3\times 3}$. Then, the above relations for turbulent energy and dissipation rate become

[TABLE]

In contrast to [21, 22], we use here a simplified energy spectrum

[TABLE]

which only reflects Kolmogorov’s 5/3-Law but ensures the differentiability of the Gaussian field and enables the corresponding distribution function $F_{E}:[\kappa_{1},\kappa_{2}]\times\mathbb{R}^{+}_{0}\rightarrow[0,1]$ and its inverse $F_{E}^{-1}:[0,1]\times\mathbb{R}^{+}_{0}\rightarrow[\kappa_{1},\kappa_{2}]$ to be stated explicitly by analytic expressions

[TABLE]

The $\xi$-dependent wave numbers $\kappa_{i}$, $(i=1,2)$, $0<\kappa_{1}<\kappa_{2}$ result from the relations in (B.1) as

[TABLE]

with $a_{E}=(54/8\xi)^{2}$, $b_{E}=(c_{E}/18^{1/3}-(128/3)^{1/3}a_{E}/c_{E})^{1/2}$, $c_{E}=(9a_{E}^{2}+\sqrt{786a_{E}^{3}+81a_{E}^{4}})^{1/3}$. The analytic expressions are advantageous in the sampling of the random field as they speed up drastically the numerical computations. Concerning the temporal correlation function we use the model of [16, 22]

[TABLE]

with $t_{L}=0.212$ being the dimensionless decay time of the turbulent structures.

Appendix C Viscous stationary uni-axial fiber jet model and numerical scheme

For stationary considerations the fiber jet model is formulated in the Eulerian (spatial) description. This is achieved by re-parametrization of the Lagrangian (material) setting of System 3 via an oriented time-dependent bijective mapping

[TABLE]

Assuming sufficient regularity, a scalar convective speed $u$ and a spatial Jacobian $j$ corresponding to $S$ is defined by

[TABLE]

In the Eulerian description the scalar speed $u$ becomes the Lagrangian multiplier to the assumption of global arc-length parametrization $\lVert\partial_{s}\mathbf{r}\rVert=1$ and is hence an additional unknown of the system.

Starting from System 3, the stationary problem in Euler description for an uni-axial fiber with orientation $\boldsymbol{\tau}=\mathbf{e_{g}}$ and purely viscous material behavior ($\mathrm{De}\rightarrow 0$) is given in non-dimensional form by System 8 on the stationary domain $\Omega=(0,1)$ with tangential aerodynamic force $f_{air}$, i.e. $f_{air}=\mathbf{f}_{air}\cdot\boldsymbol{\tau}$. The pressure of the viscous fiber jet model satisfies the relation $p=-{\sigma}/{3}$.

System 8** (Stationary viscous uni-axial fiber model).**

Kinematic equations as well as material laws in $\Omega$:

[TABLE]

Boundary conditions:

[TABLE]

System 8 is a boundary value problem of ordinary differential equations of first order on a fixed domain. For its numerical solution we employ the continuation-collocation method, which has been used successfully in [2, 3, 31]. The collocation method is a three-stage Lobatto IIIa formula, see, e.g., [15]. Mesh selection and error control are based on the residual of the continuous solution [17]. The resulting nonlinear system of equations is solved using a Newton method with numerically approximated Jacobian. This is a classical approach that is provided in the software MATLAB***For details see http://www.mathworks.com. by the routine bvp4c.m. Its applicability depends on the convergence of the Newton method that is crucially determined by the initial guess. The initial guess is adapted iteratively by means of a continuation method from [2, 3, 31]. For System 8 we introduce a single continuation parameter $c\in[0,1]$ for the viscous, gravitational and aerodynamic forces as well as for the heat exchange. In particular we set

[TABLE]

The starting solution for the continuation that belongs to $c=0$ is a stress-free straight fiber with constant speed and temperature taken from the inflow, i.e., $\sigma=0$, $T=1$ and $u=1$.