Numerical Study of Melt-Spinning Dynamic Parameters and Microstructure Development with Ongoing Crystallization

Abstract

In response to an investigation on the paths of changes in the crystallization and radial differences during the forming process of nascent fibers, in this study, we conducted numerical simulation and analyzed the changes in crystallization mechanical parameters and tensile properties through a fluid dynamics two-phase model. The model was based on the melt-spinning method focusing on melt spinning, the environment of POLYFLOW, and the method of joint simulation, coupled with Nakamura crystallization kinetics, including the development of process collaborative parameters, stretch-induced crystallization, viscoelasticity, filament cooling, gravity term, inertia, and air resistance. Finally, for nylon 6 BHS and CN9987 resin spinning, the model successfully predicted the distribution changes in temperature, velocity, strain rate tensor, birefringence, and stress tensor along the axial and radial fibers and obtained the variation pattern of fibers’ crystallinity along the entire spinning process under different stretching rates. Furthermore, we also explored the effects of spinning conditions, including inlet flow rate, winding speeds, and the extrusion temperature, on the fibers’ crystallization process and obtained the influence rules of different spinning conditions on fiber crystallization. Knowing the paths of changes in mechanical performance can provide important guidance and optimization strategies for the future industrial preparation of high-performance fibers.

1. Introduction

Melt spinning is a key technology, widely used in the industrial manufacturing of synthetic fibers. Due to the many limitations of existing technology, significant improvements in fibers leading to higher performance are constrained. In long spinning, the process of interaction between the parameters that affect the changes in the fiber structure and performance is far from adjusting a single factor that could allow for the optimization of fiber structure and performance [1,2,3]. The collaborative regulation of multiple factors brings an enormous workload to experiments, and experimental results are often limited to specific systems, with poor scalability [4,5,6]. The established online detection methods are insufficient in detecting the variation path of fibers’ microscopic parameters throughout the spinning process; however, this is crucial for understanding and predicting the formation of non-equilibrium fibers’ structures during the spinning process.

Numerical simulation is an important tool that can provide insights and predictions of the fiber structure during the formation process at the microscale, as well as motion state information, such as the velocity field, temperature field, and stress field [7,8,9]. To theoretically understand the mechanism of fiber formation during the spinning process and compensate for the shortcomings of experimental methods, Keunings et al. [10] modeled and simulated a three-dimensional spinning process of UCM fluids with the method of finite elements. The work simulated polymer melt as the isothermal Newtonian fluid in their modeling method. However, taking into account crystallization as a process of heat as well as energy transformation, the significant drawback of isothermal Newtonian fluid is that it cannot handle strain changes related to nonlinear viscoelastic effects. Furthermore, it cannot reflect the deformation of fibers during extrusion, which is inconsistent with reality. Huynh and Tanner [11] simulated and analyzed the non-isothermal spinning process of Newtonian fluids. The research results showed that there was a significant temperature difference in the radial direction of fibers, and the above simulation analysis did not consider the crystallization phase transition process of polymer materials. In fact, the equation used in the model was the Maxwell constitutive equation, which reflected the viscoelastic properties of fluids and affected the crystallization of fiber formation at length. This equation has the disadvantage of being a single-phase model and cannot reflect the formation of the neck (for example, for a winding speed of above 2000 m/min, their model’s reflected speed will not reach the platform zone). In Wilkie A’s study [12], there was a significant temperature deviation in the radial direction of cooling fibers, with the fibers’ central temperature higher than the surface temperature. Although they simulated the amorphous and semi-crystalline phases as different constitutive models, they ignored the effects of fiber cooling, air resistance, and, in some cases, inertial forces on boundary conditions, free surfaces, and phase transitions when establishing the model. Thus, the robustness of the model did not match the actual winding speed and flow rate of spinning. A similar viewpoint was also expressed in Tae Hwan Oh’s study [13], in which the fibers’ temperature was non-uniformly distributed on the cross-section.

These researchers used the average temperature of the cross-section instead of the fibers’ surface temperature, which can lead to an overestimation of thermal conductivity with the surrounding environment, resulting in incorrectly calculated temperature fields. Due to the fact that stress is a function of temperature in the flow of viscoelastic fluids, deviations in temperature field calculations can cause certain deviations in the calculated stress, leading to incorrect molecular orientation and crystallinity, ultimately affecting the evaluation of the fiber structure and properties. Sun and Wang et al. [14,15] considered the phenomenon of flow-induced the crystallization and established a simulation of two-dimensional viscoelastic spinning. The model predicted extrusion swelling, temperature, and birefringence data at the spinneret, which was consistent with experiments. However, the established model was essentially a single-phase model, which could not describe the radial structure and crystal changes of fibers. They essentially did not append different constitutive equations for the crystalline and amorphous phases, which, in a sense, limited their understanding of the complex phenomena that occur during high-melt-spinning processes. Moreover, the fibers’ solidification point is not naturally generated in the model formula but, rather, a transition point after reaching the polymer glass transition temperature by setting a fixed temperature.

In this study, based on the environment of development regarding POLYFLOW [16,17], a method of joint simulation was used to couple the Nakamura crystallization kinetics equation, and a two-phase mathematical model with a process of synergistic parameters including stretch induced the crystallization, viscoelasticity, wire cooling, gravity term, inertia, and air resistance. The model can predict the changes in crystal mechanical parameters, tensile properties, and microcrystalline crystallinity of fibers during the process of fiber formation, including radial changes in fiber temperature, velocity gradient, birefringence index, strain rate tensor, and stress tensor distribution, along the axial and radial directions, as well as the effect of different stretching on fiber crystallization. In addition, the influence of spinning conditions, including the inlet flow rate, winding speed, and extrusion temperature, on the changes in the fiber crystallization process was explored through the model. Monitoring the paths of changes in these process parameters is crucial for the industrial preparation of high-performance fibers.

2. Mathematical Model

2.1. Method

The problem of flow considered in our study is the steady-state flow of fiber materials undergoing phase transformation and crystallization during the process of melt spinning [18,19]. In the cylindrical coordinate system, it was assumed that the flow field was axisymmetric, as shown in Figure 1. The radius of the spinneret hole is r₀, the outlet speed is V₀, the temperature is T₀, the crystallinity is $\theta$, L is the distance of the filament along the process of spinning, and V_d is the stretching speed of the spinning. For the incompressible steady-state, non-isothermal viscoelastic fluids, the process of fiber-forming phase transition involves two dynamic processes: the formation and growth of crystal nuclei, which are processes of energy conversion [20,21]. At the same time, the establishment of the following models is based on incompressible viscoelastic fluids.

As shown in Figure 1, the construction of the model was based on a cylindrical coordinate system, and the origin was selected at a position with uniform radial and axial velocities. For fluids in an incompressible stable flow state, the steady-state equation of melt spinning can be solved through the three major equilibrium equations of the continuity, momentum, and energy. In the radial transmission, assuming that the axial velocity and temperature were uniform on the cross-section, and the flow was locally uniformly uniaxial extended, the expressions of the three major equations can be expressed as Equations (1)–(3).

2.1.1. Mass Balance

$$w=\rho {\displaystyle \frac{\pi {D_0}^2}{4}}{v}_z$$

(1)

where $w$ is the polymer quality, $\rho$ is the density, $D_0$ is the fiber diameter, and $V_0$ is the velocity component of the fibers along the z direction.

2.1.2. Momentum Balance

$$w{\displaystyle \frac{d{v}_z}{dz}}={\displaystyle \frac{d}{dz}}[A({\tau }_{zz}-{\tau }_{rr})]-\pi B{\mu }_0(v_z-v_d)+\rho gA+{\displaystyle \frac{\pi s}{2}}\bullet {\displaystyle \frac{d{D}_0}{dz}}$$

(2)

In Equation (2), the first term is the inertia term, the second term is the air resistance term, and B is the Bingham number, which represents the resistance of cooling air per unit length of fiber. $v_{\mathrm{d}}$ is the velocity component of the blowing air, ${\mu }_0$ is the viscosity value of the blowing air, g is the gravitational acceleration, A is the cross-sectional area of the fibers, and S is the total surface tension acting on the filament. This equation is not directly related to the constitutive equation of the amorphous and semi-crystalline phases of fibers.

2.1.3. Energy Conversion

$$\rho {\mathrm{c}}_p{v}_z{\displaystyle \frac{dT}{dz}}=({\tau }_{zz}-{\tau }_{rr}){\displaystyle \frac{d{v}_z}{dz}}-{\displaystyle \frac{4h(T-T_0)}{D_0}}+\rho \Delta h_f{v}_z{\displaystyle \frac{d{\theta }_{\mathrm{av}}}{dz}}$$

(3)

where $C_p$ is the specific heat capacity, and it is related to the spinning temperature and crystallinity. The first term in Equation (3) represents viscous dissipation, the second term represents convective heat transfer between fibers and air, h is the coefficient of thermal convective heat transfer, and the third term represents the latent heat released due to the crystallization of fibers. $\Delta h_{\mathrm{f}}$ is the crystallization heat per unit mass, and ${\theta }_{\mathrm{av}}$ is the average crystallinity of the entire spinning system. In establishing the energy equation, Denn et al. [22] neglected the radiation effect and the axial conductivity of heat while considering that the crystallization of fibers is an exceptional process. The last term in the energy equation retained a term of heat released by late heat.

2.1.4. Tensile Kinematics and Stress Field

Assuming that the shear stress generated in the spinneret was already in a relaxed state, and the flow field during the process of spinning was a locally uniform uniaxial tensile flow, relevant researchers have adopted this assumption during melt spinning [23,24,25]. Therefore, under the condition of constant density, the velocity tensor can be expressed as in Equation (4).

$$\nabla v=\left[\begin{array}{ccc}{\displaystyle \frac{\mathrm{d}{v}_z}{dz}}& 0& 0\\ 0& -{\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathrm{d}{v}_z}{dz}}& 0\\ 0& 0& -{\displaystyle \frac{1}{2}}{\displaystyle \frac{\mathrm{d}{v}_z}{dz}}\end{array}\right]$$

(4)

The stress tensor can be written as Equation (5).

$$\tau =\left[\begin{array}{ccc}{\tau }_{\mathrm{zz}}& & \\ & {\tau }_{rr}& \\ & & {\tau }_{\theta \theta }\end{array}\right]$$

(5)

where $zz$, rr and $\theta \theta$ represent stress components in three different directions.

2.1.5. Constitutive Model

Based on the process of spinning design and the characteristics of the materials, the characteristics of the constitutive equations of the Maxwell model, Giesekus model, and Phan Thien Tanner (PTT) were compared. The expression of control equations for the three constitutive models is (6–18).

Maxwell model:

$$\sigma +\lambda [{\displaystyle \frac{\partial \sigma }{\partial t}}+v\bullet \nabla \sigma -(\nabla v)\bullet \sigma -\sigma \bullet {(\nabla v)}^T-\sigma \bullet (v\nabla \ln T)]=\mu [\nabla v+{(\nabla v)}^T]$$

(6)

where $\sigma$ is the stress tensor, $\mu$ is the shear viscosity, t is the melting time, and $\lambda$ is the relaxation time.

For the flow of incompressible steady-state fluid, the above equations can be written as a component form with the following expressions (7) and (8) during simulation [13].

$${\tau }_{\mathrm{zz}}+\lambda [v_z{\displaystyle \frac{d{\tau }_{zz}}{dz}}-2{\tau }_{zz}{\displaystyle \frac{d{v}_z}{dz}}-v_z{\tau }_{zz}{\displaystyle \frac{d\ln T}{dz}}]=2\mu {\displaystyle \frac{d{v}_z}{dz}}$$

(7)

$${\tau }_{\mathrm{rr}}+\lambda [v_z{\displaystyle \frac{d{\tau }_{rr}}{dz}}-2{\tau }_{rr}{\displaystyle \frac{d{v}_z}{dz}}-v_z{\tau }_{rr}{\displaystyle \frac{d\ln T}{dz}}]=2\mu {\displaystyle \frac{d{v}_z}{dz}}$$

(8)

where ${\tau }_{\mathrm{zz}}$ is the axial force of the fiber, and ${\tau }_{\mathrm{rr}}$ is the radial force of the fiber.

The mathematical relationship among $\sigma$ and ${\tau }_{\mathrm{zz}}$, ${\tau }_{\mathrm{rr}}$ is Equation (9).

$$\sigma ={\tau }_{\mathrm{zz}}-{\tau }_{rr}$$

(9)

Giesekus model:

$$\tau ={\displaystyle \sum _{\mathrm{t}}^n{\tau }_j}$$

(10)

$${\displaystyle \frac{\partial {\tau }_{\mathrm{j}}}{\partial t}}+v\bullet \nabla {\tau }_j-\nabla v^T\bullet {\tau }_j-{\tau }_j\bullet \nabla v+{\displaystyle \frac{1}{{\lambda }_j}}({\tau }_j+{\displaystyle \frac{\psi {\lambda }_j}{{\eta }_j}}{\tau }_j)={\displaystyle \frac{{\lambda }_j(\nabla v+\nabla v^T)}{{\eta }_j}}$$

(11)

where ${\lambda }_j$ and ${\eta }_j$ are the relaxation time and shear viscosity under the j model; $\psi$ is the liquidity factor, with a value range of 0–1.

PTT model:

$$\exp [{\displaystyle \frac{\epsilon \lambda }{\eta }}\mathrm{tr}(\tau )]\tau +\lambda [(1-{\displaystyle \frac{\delta }{2}})\stackrel{\nabla }{\tau }+{\displaystyle \frac{\delta }{2}}\stackrel{\Delta }{\tau }]=2\eta D$$

(12)

wherein

$$\stackrel{\nabla }{\tau }=v\bullet \nabla \tau -\nabla v^T\bullet \tau -\tau \bullet \nabla v$$

(13)

and

$$\stackrel{\Delta }{\tau }=v\bullet \nabla \tau +\tau \bullet \nabla v^T+\nabla v\bullet z$$

(14)

The above equations, Equations (13) and (14), represent the upper- and lower-body derivatives of the polymer stress tensor $\tau$ under the condition of a steady state, respectively. $\delta$ is the slip factor between the continuous medium and the molecular network model, $\epsilon$ represents the change in tensile viscosity, $\lambda$ is the relaxation time, and D is the strain rate tensor; its expression is shown in Equation (15).

$$D={\displaystyle \frac{(\nabla v+\nabla v^T)}{2}}$$

(15)

Assuming that the principle of “time temperature equivalence” conformed to the temperature-dependent relationship of the viscoelasticity of spinning materials, the relationship between viscosity $\eta$, relaxation time $\lambda$, and temperature T can be expressed as (16) and (17).

$$\eta (T)={\alpha }_T(T_0)$$

(16)

$$\lambda (T)={\alpha }_T(T_0)$$

(17)

where $T_0$ is the reference temperature, T is the absolute temperature, and ${\alpha }_T$ is the conversion factor; its value can be calculated using the Arrhenius equation, Equation (18).

$$\ln ({\alpha }_T)=\ln ({\displaystyle \frac{\eta }{{\eta }_0}})={\displaystyle \frac{\Delta H}{R}}({\displaystyle \frac{1}{T}}-{\displaystyle \frac{1}{T_0}})$$

(18)

where $\Delta H$ is the viscous activation energy of the material, and R is the gas constant at the reference temperature.

2.1.6. Definition of Dimensionless Variables

To improve computational efficiency, the dimensionless variable expressions (19–31) were introduced.

Inertia:

$$D_1={\displaystyle \frac{\rho {v_0}^2}{G}}$$

(19)

Air drag:

$$D_2={\displaystyle \frac{\mathsf{\pi }\eta LB\rho {v_0}^2}{GW}}$$

(20)

Gravity:

$$D_3={\displaystyle \frac{gL\rho }{G}}$$

(21)

Surface tension:

$$D_4={\displaystyle \frac{S}{2G}}\sqrt{{\displaystyle \frac{\mathsf{\pi }\rho v_0}{W}}}$$

(22)

Heat convection:

$$D_5={\displaystyle \frac{2Lh}{c_p}}\sqrt{{\displaystyle \frac{\mathsf{\pi }}{\rho v_{0}W}}}$$

(23)

Viscous dissipation:

$$D_6={\displaystyle \frac{G}{\rho {\mathrm{c}}_p{T}_0}}$$

(24)

Latent heat of crystallization:

$$D_7={\displaystyle \frac{\Delta H_{\mathrm{f}}{\theta }_{\infty }}{c_p{T}_0}}$$

(25)

Dimensionless axial distance:

$${\mathrm{z}}^{*}={\displaystyle \frac{z}{L}}$$

(26)

Dimensionless velocity:

$${v_z}^{*}={\displaystyle \frac{v_z}{v_0}}$$

(27)

Dimensionless temperature:

$$T^{*}={\displaystyle \frac{T}{T_0}}$$

(28)

Dimensionless extra stress tensor:

$${\tau }^{*}={\displaystyle \frac{\tau }{G}}$$

(29)

Relative velocity:

$$v_r={\displaystyle \frac{v_d}{v_0}}$$

(30)

Deborah number of semi-crystalline phase:

$$D_{\mathrm{esl}}={\displaystyle \frac{\lambda \bullet v_0}{L}}$$

(31)

2.1.7. Evolution Methods

There are two types of variables in the model; one includes the macroscopic velocity change V, the temperature change T, and the fiber cross-sectional size change. Another type is the microstructural changes in total stress tensor S, strain rate tensor D, and birefringence $\Delta n$ in the amorphous and semi-crystalline regions. Before beginning crystallization, the independent fluid variable is V as well as T. D_zz is relative to the axial direction of the fibers. D_rr is relative to the radial direction of the fibers.

Momentum balance:

Both sides of Equation (2) use $L/W{v}_0$, which can eliminate the diameter $D_0$ of the filament, and, at this point, the momentum balance control equation will become Equation (32).

$$D_1{\displaystyle \frac{\mathrm{d}{v_z}^{*}}{{d_z}^{*}}}={\displaystyle \frac{d}{d{z}^{*}}}[{\displaystyle \frac{{{\tau }_{zz}}^{*}-{{\tau }_{rr}}^{*}}{{v_z}^{*}}}]-D_2({v_z}^{*}-v_r)+{\displaystyle \frac{D_3}{{v_z}^{*}}}-D_4{({v_z}^{*})}^{-{\displaystyle \frac{3}{2}}}{\displaystyle \frac{d{v_z}^{*}}{d{z}^{*}}}$$

(32)

Energy equation:

Multiply both sides of Equation (3) by $L{T}_0/\rho c_p{v}_0$; thus, after sorting, Equation (33) can be obtained.

$${\displaystyle \frac{\mathrm{d}{T}^{*}}{d{z}^{*}}}=-D_5(v_z^{*}{)}^{-{\displaystyle \frac{1}{2}}}(T^{*}-T_0)+D_6{\displaystyle \frac{{{\tau }_{zz}}^{*}-{{\tau }_{rr}}^{*}}{{v_z}^{*}}}{\displaystyle \frac{d{v_z}^{*}}{d{z}^{*}}}+D_7{\displaystyle \frac{1}{{v_z}^{*}}}{[-\ln (1-x)]}^{{\displaystyle \frac{m-1}{m}}}(1-x)\exp (\mathrm{tr} {\tau }^{*})$$

(33)

The evolution equation for birefringence is as follows:

$${\displaystyle \frac{\mathrm{d}\Delta n^{*}}{d{z}^{*}}}=(1-\theta ){\displaystyle \frac{c_{op}\mathsf{\mu }{v}_0}{L_0}}({{\tau }_{zz}}^{*}-{{\tau }_{rr}}^{*})+\theta f_c{\displaystyle \frac{\mathrm{d}\Delta {n_c}^{*}}{d{z}^{*}}}$$

(34)

where $f_c$ is the crystal orientation factor, $f_c=0.9$ [16], and $c_{op}$ is the optical coefficient of stress.

The evolution equation for $S_{zz}$ is as follows:

$${\displaystyle \frac{{\mathrm{dS}}_{\mathrm{zz}}}{d{z}^{*}}}={\displaystyle \frac{2S_{zz}}{{v_z}^{*}}}{\displaystyle \frac{d{v_z}^{*}}{d{z}^{*}}}-{\displaystyle \frac{\sigma S_{zz}}{{v_z}^{*}{D}_{esc}}}+{\displaystyle \frac{2}{3v_{z*}}}{\displaystyle \frac{d{v_z}^{*}}{d{z}^{*}}}-{\displaystyle \frac{1.7(1-w)}{{v_z}^{*}}}({\displaystyle \frac{7}{15}}+{\displaystyle \frac{11}{14}}{S}_{zz}){\displaystyle \frac{d{v_z}^{*}}{d{z}^{*}}}-{\displaystyle \frac{3w{S}_{zz}}{{v_z}^{*}}}(S_{zz}+0.34){\displaystyle \frac{{\mathit{dv}}_z^{*}}{d{z}^{*}}}$$

(35)

The $S_{\mathrm{rr}}$ and $S_{\theta \theta }$ components are eliminated using axial symmetry and the condition for S being traceless.

The expression before crystallization starts is (36).

$${\tau }^{*}=E{\mathrm{c}}^{*}-\delta$$

(36)

The expression after crystallization starts is (37).

$${\tau }^{*}={\displaystyle \frac{E}{1-x}}{c}^{*}-\delta +3S+6D_{esc}{({\nabla }^{*}{v}^{*})}^T:<SSSS>$$

(37)

Considering the influence of the elastic modulus of the amorphous phase on the crystalline phase, the influence factor $E/(1-x)$ is introduced in Equation (37).

Equation (38) can be obtained through the energy equation.

$$\rho c_p({\displaystyle \frac{\partial T}{\partial t}}+v_z{\displaystyle \frac{\partial T}{\partial t}})={\displaystyle \frac{1}{r}}{\displaystyle \frac{\partial }{\partial r}}(r{\lambda }_{*}{\displaystyle \frac{\partial T}{\partial r}})+{\eta }_0{r}^2+Q$$

(38)

where $\rho$ is the density, $C\mathrm{p}$ is the specific heat capacity at time t, t is the time, ${\eta }_0$ is the kinematic viscosity, T is the spinning temperature, r is the axial radius, z is the axial length, ${\lambda }_{*}$ is the thermal conductivity, $v_z$ is the flow velocity along the z direction, and $Q^{*}$ is the release of heat per unit volume during the crystallization process, which is the latent heat term of crystallization.

Among them, the latent heat of the crystallization term $Q^{*}$ can be expressed by the crystallization kinetics expression (39) [26].

$$Q^{*}=\rho {\theta }_{\infty }\Delta H_C{\displaystyle \frac{\mathrm{d}\theta }{dt}}$$

(39)

where ${\theta }_{\infty }$ is the maximum crystallinity, and $\Delta H_c$ is the enthalpy of crystallization. It is not difficult to see that crystallization is not only affected by the extrusion temperature, cooling air temperature, and tensile stress but that these variables also change over time.

Norioki K et al. [27], H W Jing et al. [28], and Evans et al. [29] each derived the dependence of crystallinity $\theta$ on time from a statistical perspective, and their expression equations are shown in (40).

$${\displaystyle \frac{\mathrm{d}\theta }{dt}}=n{k}_s(T,E(t))(1-\theta ){[-\ln (1-\theta )]}^{{\displaystyle \frac{n-1}{n}}}$$

(40)

Among them, n is the Avrami index, which is used to describe the mechanism and rate of crystal growth during the process of crystallization. Its value can be fitted by the Avrami equation experimental data. When the value of n is 1, it indicates that the crystal is growing linearly in one dimension. When the value of n is 2, it indicates that the crystal is growing in two dimensions. When the crystal is growing in three dimensions, the value of n is 3.

By using the iterative integration method on Equation (40) above, an integrated equation can be obtained, as shown in Equation (41).

$$\theta =1-\exp [-({\displaystyle \underset{0}{\overset{t}{\int }}{k}_s(T,E(t)\mathrm{dt})}{)}^n]$$

(41)

where $k_s(T,E(t))$ is the crystallization rate at time t.

For flow (shear)-induced crystallization, the rate of crystallization is controlled by Equation (42).

$$k_s(T,\stackrel{\bullet }{\gamma }){=\left(\ln 2\right)}^{{\displaystyle \frac{1}{n}}}({\displaystyle \frac{1}{{\kappa }_t}})\exp ({\displaystyle \frac{-\Delta H}{R(T-T_{\infty })}})\bullet \exp (-{\displaystyle \frac{k_k}{Tf(T_m(T,\stackrel{\bullet }{\gamma })-T)}})$$

(42)

where $T_{\infty }$ represents the temperature at which the filament glass melts along the spinning process, ${\kappa }_{\mathrm{t}}$ is the reference factor that affects crystallization, $\stackrel{\bullet }{\gamma }$ is the shear rate, $k_k$ is the nucleation index, $T_{\mathrm{m}}$ is the melting point, and R is the gas constant.

The control equation for the crystallization coupled with Nakamura kinetics [30] induced by stretching (orientation) is (43).

$$k_s(T,f_{ac})=k_{\max }\exp [-4\ln 2{\displaystyle \frac{T-T_{\max }}{D^{*}}}+{\mathit{cf}}_{ac}^2]$$

(43)

where $f_{ac}$ is the orientation factor, c is the stress-inducing factor, $D^{*}$ is the half-crystalline width, $K_{\max }$ is the maximum crystallization rate, and $T_{\max }$ is the temperature at the maximum rate of crystallization.

Almost all crystallization processes in chemical fiber technology are carried out in an oriented state [31]. Therefore, the rate of the crystallization kinetics model constructed in our models is based on the above equations.

• wherein

$$f_{ac}={\displaystyle \frac{\Delta \alpha }{\Delta {\alpha }_i}}$$

(44)

$\Delta \alpha$ is the birefringence in the amorphous region, and the calculation equation is (44). $\Delta {\alpha }_{\mathrm{i}}$ represents the inherent birefringence of the corresponding amorphous region. The amorphous birefringence is related to the laws of stress optics and the first normal force difference.

$$\Delta \alpha =c_{op}^{*}\eta {\displaystyle \frac{v_0}{L_0}}({\tau }_{zz}-{\tau }_{rr})$$

(45)

where $c_{op}^{*}$ is the stress optical coefficient, and ${\tau }_{zz}$ and ${\tau }_{\mathrm{rr}}$ represent the stress components along the z and r directions, respectively. The above three equations are the final kinetic models for evaluating crystallization.

Viscosity function:

For the flow of an incompressible polymer, the expression for its viscosity function can be obtained using the method of Sun et al. [14] as (46).

$$\eta (T)=T\bullet \exp \langle {\displaystyle \frac{\Delta H}{R}}({\displaystyle \frac{1}{T}}-1)\rangle$$

(46)

In the formula, R is the gas constant at the standard atmospheric pressure. Considering the influence of crystallinity on viscosity, the dependence of viscosity on crystallinity is further quantified by introducing Equation (47).

$${\eta }^{*}(T)=\eta (T)\exp \langle 4{\theta }^2\rangle$$

(47)

2.2. Control of Boundaries

${\Gamma }_{\mathrm{i}n},{\Gamma }_{wall},{\Gamma }_{out}$, and ${\Gamma }_{free}$ represent the inlet boundary of the fluid, the wall of the spinneret, the outlet of the fluid, and the free surface of the fibers. The boundary condition for the spinning process was defined as (48–51), taking into account the full development of the fluid motion.

$${\Gamma }_{\mathrm{i}n}: v=v(r),{\tau }_p={\tau }_p(r)$$

(48)

$${\Gamma }_{\mathrm{wall}}: v=0,T=T_{wall}$$

(49)

$${\Gamma }_{\mathrm{out}}: v_z=v,F_r=0,\nabla T\bullet n=0$$

(50)

$${\Gamma }_{\mathrm{f}ree}: \sigma \bullet n=0,v\bullet n=0,\nabla T\bullet n=-Bi(T-T_0)$$

(51)

where $\sigma$ is the total stress tensor and Bi is the Biot number; its calculation equation is shown in (52).

$$B\mathrm{i}\equiv {\displaystyle \frac{f_h\bullet L}{k}}$$

(52)

where $f_h$ is the thermal conversion efficiency, calculated using the method of Kase and Matsuo [32] and expressed as (53).

$$f_h=1.98A^{-0.333}{v}^{0.333}{[1+{\displaystyle \frac{8v_{\mathrm{air}}}{v}}]}^{0.167}$$

(53)

A is the cross-sectional area of the free surface of the fibers, and $v_{\mathrm{air}}$ is the speed of blowing air.

3. Numerical Method and Solution Procedure

During the simulation, the commercial software POLYFLOW was specifically used for solving non-Newtonian fluids. The model of crystallization kinetics was defined based on the transport characteristics of POLYFLOW and user-defined functions. The steps of the specific solution are shown in Figure 2. The hardware configuration environment used Windows 10 as the operating system, with an Inter (R) Core (TM) i9-10900, a CPU frequency of 208 GHz, and 32 GB of RAM on a PC. Considering the expansion of the fiber outlet and ensuring the accuracy of the calculation, mesh refinement was carried out at the inlet of fluid, followed by a gradually sparse mesh with a sparsity deviation of 2.05. In each unit, the velocity vector and stress component were solved using the DEVSS SU difference function [19].

The problem of convergence is the key to the success of simulation calculations. In our study, in order to improve the possibility of convergence and ensure the accuracy of simulation calculations, iterative functions and linear functions were used to solve the stress and strain rate tensors in the control equations. Considering that the deviatoric stress tensor is decomposed into elastic and viscous components using the EVSS (Elastic Viscous Splitting of the Stress) method, the SUPG method was advantageous for coupling between equations, with a numerical convergence standard of 10⁻⁶.

4. Numerical Simulations and Input Parameters

Selecting circular cross-section fibers as simulation objects, using Gambit to establish a geometric model of circular cross-section fibers, and meshing the geometric model, as shown in Figure 3, all meshes in the model geometry were tetrahedral structures. For the fluid inlet, surface mesh generation was used, and the free end was swept. To ensure the accuracy of simulation calculations, the “Implicit Euler” method was coupled with the “Picard Iterative” function. For the coupling calculation between free surface and fluid equations, in order to avoid the influence of free boundary repositioning on some grids, the “Upwinding in the kinematic equation” and “Optotimesh-3D” grid reset methods were adopted. The maximum ratio of aspect for grid partitioning was 3.328, the minimum unit mass was 0.021, the minimum orthogonal mass was 0.11, the maximum deviation was 0.102, and the shrinkage difference was 0.002. The total number of vertices in the meshing situation was 1035, the number of data segments was 1858, and the total number of planes was 824.

During the cooling process of fiber formation, when the surface temperature of the fibers is lower than the melting point temperature, the fibers will undergo crystallization, and different conditions of spinning and winding speeds will cause radial differences in the fibers and changes in the path of microcrystalline crystallization. Therefore, we simulated and analyzed mathematical models under the process of radius r₀, spinning process L, and temperature T at the outlet of the spinneret.

4.1. The Influence of Different Constitutive Models on Extrusion

Three models, the Maxwell model, PTT model, and Giesekus model, were selected based on experience and the theory of fluid [33,34]. The Maxwell model uses two parameters, shear viscosity and relaxation, to represent the liner viscoelastic behavior of fluids. However, this model cannot accurately describe the viscoelastic properties of fluids during deformation and flow. The PTT model was the constitutive equation derived from network theory. What’s more, it can better describe the rheological properties of polymer melts as well as high shear rates and deformation tensors. The Giesekus model can describe the region where the viscosity of a fluid remains constant at low and high shear rates. The forming process of fibers is mainly influenced by stretching and cooling air, and the crystallization process is actually a combination of static crystallization and stretch-induced crystallization. The process of non-isothermal crystallization is essentially a process of energy conversion. Therefore, the path change of temperature parameters plays a crucial role in the crystallization of fibers. The energy equation was simulated using the material dataset publicly available in Nylon 6. At the same time, there was a maximum temperature difference between the core of the spinneret and the surface, which gradually decreased throughout the process of spinning. The variation throughout the process of spinning is shown in Figure 4.

Figure 4 indicates the predicted temperature distribution along the fiber centerline and free surface. Due to the effect of thermal convection cooling, the surface of fibers and air undergo convective heat transfer. When fluid particles move downstream, the temperature decreases, and the temperature along the surface of fibers is lower than that along the center of fibers. Considering that viscous heat generation is ignored in the calculation, combined with the data values in Figure 4a. The temperature predicted by the PTT model differs from the maximum experimentally measured value by less than 3%. Therefore, PTT shows the best agreement among the three models. At the same time, the PTT model is the constitutive equation that describes the viscoelasticity of the polymer using molecular network models and relaxation time. The model contains two significant parameters, that is, the stretch viscosity control parameter as well as the continuous medium slip control parameter.

By analyzing Figure 4b, it can be observed that there is a significant temperature difference between the center of fibers and the outer surface a certain distance from the extrusion nozzle. Moreover, the temperature difference between the surface of fibers and the center of fibers at a speed of 15 m/s was greater than that at 70 m/s, indicating that the lower the stretching rate, the greater the difference between the fiber surface temperature and the center of temperature.

4.2. Input Parameters

In our study, based on the established kinetic model, Nylon 6 BHS and CN9987 resin materials were used as the carrier.

Table 1. Geometry and model parameters for dynamic equations.

Processing Parameter	Value
Diameter of entrance channel, D0 (mm)	0.76
Lenth of exit channel, L0 (mm)	2.42
Distance from spinneret, L (mm)	2 × 103
Extrusion temperature, T0 (k)	533
Molten temperature, Tm (k)	503
Maximum rate temperature, Tmax (k)	419
Cooling temperature, Tc (k)	298
Blowing speed, Vair (m/s)	0.4

,

Table 2. Physical and rheological parameters.

Properties	Parameters	References
Viscosity, $\eta$ (Pa·s)	750	[14]
Relaxation time, $\lambda$ (s)	0.015	[31]
Elongational parameter, $\epsilon$	0.2	[31]
Mobility factor, $\xi$	0.001	[31]
Additional viscosity, (Pa·s)	0.25	[31]
Density, $\rho$ (kg/m3)	973	[15]
Thermal conductivity, k (w/m/k)	0.2	[31]
Heat capacity, C (J/kg/k)	1.3497 × T + 1982.76	[31]
Arrhenius, $\alpha$ (k)	4470	[31]
Ultimate crystallinity, ${\theta }_{\infty }$	0.45	[31]
Maximum crystallization rate at Tmax, Kmax (s−1)	0.14	[31]
Crystallization half width, $D_{\left(\mathrm{k}\right)}^{*}$	46	[15]
Stress induced coefficient, c	225	[15]
Stress-optical coefficient, $C_{\mathrm{op}}^{*}$ (m2N−1)	1.25 × 10−9	[15]
Amorphous intrinsic birefringence, $\Delta {\alpha }_{\mathrm{i}}$	8.25 × 10−2	[31]
Heat of fusion, Hf (Jkg−1)	1.884 × 105	[15]

and

Table 3. Conditions of melt spinning in the simulation.

Number	Volume Flow Rate (mm3/s Hole)	Spinning Temperature (k)	Circular Air Temperature (k)	Spinning Speed (m/min)
1	30	533	293	3000
2	45	533	293	3000
3	60	533	293	3000
4	75	533	293	3000
5	60	533	293	1800
6	60	533	293	2400
7	60	533	293	3600
8	60	533	293	4200
9	60	535	293	3000
10	60	538	293	3000
11	60	542	293	3000
12	60	545	293	3000
13	60	533	298	3000
14	60	533	300	3000
15	60	533	304	3000

, respectively, list the physical as well as rheological parameters, process parameters, and operational parameters of the simulated input materials, corresponding to the parameters in the simulation. These parameters are derived from data reported by researchers [13,15], so these data were cited in our study.

5. Results and Discussion

The uneven distribution of melt-spinning fibers in the radial structure is a key factor restricting the development of high-performance, high-value-added fibers. Throughout the entire process of melt spinning, the differences in the radial structure distribution of fibers and the changes in paths of crystallization are the coordinated effects of many parameters. It is difficult to achieve effective control through unilateral parameter control, and it is not possible to discover the regular changes and quantitative value changes in the path. Even existing online detection methods such as infrared detection and high-speed photography can only perform discrete detection on a limited number of fixed positions under the condition of characteristic environments, and the results of detection are not only affected by subjective observations but also greatly interfered by the accuracy of detection technology.

Therefore, the model we established above was used to simulate the process of melt-spinning fiber formation, taking crystallization into account. The joint simulation method was used to predict the temperature, velocity, strain rate tensor, and stress tensor distribution of nylon 6 fibers along different directions. Meanwhile, through the model, we also explored the effects of conditions of spinning, including the inlet flow rate Q, winding speed V, extrusion temperature T, and crystallization of fibers $\theta$. The model is more applicable to “Melt-spinning Methods” and “Viscoelastic Non-Newtonian Fluids”. The specific discussion is as follows.

5.1. Research on the Paths of Changes in Speed

The new model broke through the traditional assumption of radial temperature and uniform molecular orientation; introduced free surface thermal convection B, and semi-crystalline width $D_{\left(\mathrm{k}\right)}^{*}$; and established crystalline and semi-crystalline phases as the Giesekus and PTT constitutive models. At the same time, the Nakamura dynamic equation was introduced in the crystallization model, so the crystallization inflection point was that the appearance was no longer influenced by human intervention but generated in a natural state through established models.

To verify the superiority of our model, we compared the predicted speed results of our model with the traditional single-phase model mentioned in the Introduction [15]. Under the same process of parameters, we compared the numerical simulation results of the stretching speed distribution of our model at a winding speed of 30 m/s with the experimental results [37], which were consistent. This indicated that our method simulated the forming process of fibers in a similar manner to the actual spinning and the process of fiber production, providing credibility for the model parameters used in further calculations. The results are shown in Figure 5.

Figure 6 indicates that five speeds represent the changes in high-speed spinning, and it was found that the trend of change in the stretching rate is the same at different stretching rates. However, the magnitude of the change in the stretching rate was variable. When L =1 m, the distribution of stretching rate values is 23.7 m/s, 29.9 m/s, 35.8 m/s, 41.4 m/s, and 46.8 m/s. The faster the winding speed, the more pronounced the S-shaped distribution of the stretching rate of the melt, and there was a region of rapid increase, followed by a stable speed. At the same time, it can be seen from Figure 6b that at a distance of 0.1 mm from the spinneret, the radial velocity component along the z-direction shows a sharp downward trend. However, as the distance from the spinneret reaches 1 mm, the velocity is relatively calm and shows a slow upward trend overall.

5.2. Temperature Distribution along the Radial Direction at Different Z Points of the Filament

The temperature distribution along the radial direction at different z points under the parameters of the wall temperature T = 533 k, stretching rate V = 30 m/s, cooling temperature T₀ = 293 K, blowing speed V₀ = 0.4 m/s, and inlet flow rate of 4.5 × 10⁻⁸ m³/s is shown in Figure 7.

Figure 7 indicates the radial temperature distribution at different axial positions in the initial deformation zone and near the solidification zone. As a result of the large temperature difference, the radial temperature distribution bends in the initial deformation zone and becomes gentle near the solidification point. At z = 5 mm, the temperature is still close to the extrusion temperature $T_{\mathrm{o}}\left(\mathrm{k}\right)$. As z increases, due to heat transfer between the surrounding medium and the surface of fibers, the temperature decreases and $\partial T/\partial \mathrm{r}$ increases. At z = 1304 mm, the leading order of radial temperature variation developed, and $(\partial T/\partial \mathrm{r})/(\partial T/\partial z)$ is 0. Further downstream, for z > 1304 mm, this leading order radial variation, $\partial T/\partial \mathrm{r}\gg \partial T/\partial z$, is shown in Figure 7b.

5.3. The Effects of the Draw-Down Velocities

To investigate the effect of different stretching speeds on the crystallization of microcrystals inside fibers, different winding speeds were applied at the end of the fibers. The results of simulation are shown in Figure 8. As shown in Figure 8, the predicted results of the model are consistent with the experimental measurements throughout the entire spinning process. Near the spinneret hole, significant fiber swelling and necking are observed due to sudden changes in boundary conditions and viscoelastic fluid properties. Under the studied conditions, the fiber diameter reached its maximum value of R = 1.032 mm at a distance of z = 0.7 mm from the spray hole on the free surface, with a swelling of 0.27 mm. And as the speed increases, the swelling gradually decreases, because the higher the speed, the greater the tensile stress on the fibers.

Figure 9 shows the distribution of the stress tensor and phase transition rate tensor during high-speed spinning in the process of fiber crystallization. As shown in Figure 9a, with an increase in the drafting speed, the average stress on the fibers increases, and the change in stress along the z direction of the spinning process shows an s-shaped trend. T₀ investigate the radial variation in total additional stress in fibers at a stretching speed of v = 30 m/s, as shown in Figure 9b, it can be observed that at z = 30 mm, the total additional tensor almost followed a straight line along the radial direction, indicating minor variation. However, as the process of spinning increases, at z = 70 mm, the total additional tensor shows an almost upward trend along the radial direction, and Figure 9c indicates the variation in stress tensor components along the z direction. Figure 9d shows the variation in the fiber strain rate tensor. As the stretching speed increases, the deformation rate in the neck area becomes sharper with an increase in flow velocity.

Figure 10 indicates the average birefringence distribution of the predicted fiber cross-sections at different z points throughout the process of spinning using numerical solutions. At a given stretching rate, the birefringence index increases with an increase in the spinning distance z and tends to reach equilibrium further along the spinning distance. For a given cross-sectional position z, the birefringence index increases with an increase in the stretching speed.

What’s more, according to Figure 10, the birefringence experimental data at three different stretching rates of 1800 m/min, 2400 m/min, and 3600 m/min were compared with the values of data simulated by the model established in this study. Firstly, from a qualitative analysis perspective, its overall trend of change remains consistent. Secondly, from a quantitative perspective, the refractive index values of fibers undergoing the same process of spinning are basically consistent. However, the consistency between the predicted values and experimental observations is not as good quantitatively as the temperature changes on the free surface, especially at high stretching rates. The reason for this is due to the lack of elongation viscosity data, and the PTT model parameters are determined by fitting shear viscosity and combining it with shear control parameters. However, the results of previous numerical studies indicate that the predicted fiber surface matches well with the experimental observations, which provides credibility for the calculation of the model. Meanwhile, the speed variation throughout the process of spinning is highly consistent with the experimental results. Combining Equations (41) and (44), birefringence is a function of temperature and stretching, and crystallization is a function of birefringence. This indirectly proves the correctness and accuracy of the model of crystallization kinetics established in this study.

Figure 11a shows the average crystallinity of the fiber cross-section throughout the process of spinning at different drafting speeds. Under the given process conditions, an increase in crystallinity with an increase in the process of spinning z ultimately reached equilibrium. The position of the equilibrium point is related to the magnitude of the stretching speed. With an increase in the stretching speed, the relative residence time of macromolecules inside the fibers will decrease, and the crystal nucleus will not have enough growth time, causing crystallization to migrate from upstream to downstream. Figure 11b shows the variation pattern of stretching speed and final crystallinity. Under the same process parameters, the higher the stretching speed, the greater the tensile stress on the fibers, and the crystallization rate is higher. On the contrary, due to the influence of tensile force, an increase in stretching speed will result in a decrease in the final fiber crystallinity. Figure 11c shows the effect of different spinning flow rates on the process of fiber crystallization. It can be observed that as the inlet flow rate increases, the crystallization rate of fibers slows down in the initial stage, but the final crystallinity of fibers increases. Figure 11d shows the effect of wall temperature on the crystallization of fibers. As the wall-melting temperature increases, the crystallization rate actually decreases, and the crystallinity of the resulting fibers also decreases.

Figure 12 summarizes the effect of $X/X_{\max }$ on the radial distribution of fibers. X represents temperature T, stress tensor $\tau$, crystallinity $\theta$, and birefringence $\Delta n$. To define ${(X/X_{\max })}_{\min }$ as the minimum value of $X/X_{\max }$, $\Delta X$ can be regarded as the quantitative analysis value for evaluating the non-uniformity of X, and the calculation expression is shown in (54).

$$\Delta X=1-({\displaystyle \frac{X}{X_{\max }}}{)}_{\min } \mathrm{at} z=L$$

(54)

If the value of $\Delta X$ is 0, it represents a uniform distribution, and if the value of $\Delta X$ is 1, it represents 100% relative variation. Regardless of the temperature, stress tensor, birefringence index, and crystallinity, the changes $X/X_{\max }$ intersect at a point of radial z = 0.016 mm. This intersection indicates that the fluid non-uniformity of the fibers during formation is minimized at this point, and the system of formation reaches the most stable state. Before 0.016 mm, the temperature T exhibits the smallest non-uniformity change. After the intersection point, these non-uniformities increase with an increase in the stretching rate, which has a significant impact on the fiber microstructure. This requires special consideration in fiber preparation and the process of control.

6. Conclusions

Through the two-phase model established in this research, which included changes in the fiber radial structure and microcrystalline crystallization in melt spinning, the change paths of the process of parameters related to the entire spinning process and fiber mechanical properties can be obtained. Compared with traditional models, firstly, the innovation of the model established in this study is the specific distinction between amorphous and semi-crystalline phases in the rheological response of the system. In addition, another key innovation of the model is that the emergence of crystallization inflection points is no longer a result of human intervention but, rather, occurs in a natural state, which is not present in traditional models. Finally, our model can accurately predict the changes in variable parameters related to crystallization, such as stress tensor, strain rate tensor, and crystallization rate, under different stretching rates. Understanding the path changes of these parameters can provide important guidance and optimization for the future industrial preparation of high-performance fibers.

The model was modified and validated through nylon 6 BHS resin and CN9987 resin spinning simulation and experimental results, and the final results showed that the model exhibited strong robustness across a wide range of spinning conditions and parameter processes. The model successfully predicted the axial and radial variations in the wire temperature, velocity, strain rate tensor, birefringence index, and stress tensor under different drafting conditions. What’s more, the influence of spinning conditions on the fiber crystallization process was studied through model simulation. The results indicated that along with an increase in the inlet flow rate, the crystallinity of the fiber increased, while with an increase in the winding speed, the crystallinity of the fiber decreased. At the same time, the wall of the spinneret hole also had an impact on the fiber crystallization, with an increase in the temperature and a decrease in the final fiber crystallinity.