Robust Numerical Resolution of Nakamura Crystallization Kinetics

The numerical prediction of crystallization transformation is of great interest in several applications. One such application is the polymer-forming process. In this short communication, the integration of the widely used Nakamura kinetics is discussed. A robust time integration method is proposed. In order to overcome its singularities, the Nakamura function is thresholded. A convergence analysis provides guidelines for the threshold values and time discretization.


Introduction
In recent decades, the numerical simulation of manufacturing processes of polymers and polymer matrix composites has been significantly improved. The latest simulation softwares [1,2,3] can account for several complex coupled non-linear physical phenomena [4]. With the development of these predictive tools, the industry could better understand, control, and eventually optimize the forming processes. today, optimizing the processing windows, in terms of process parameters or material properties, requires quantitative simulation tools. Thus an accurate modeling of fine physical phenomena is needed.
Industrial thermoplastic polymers are often semi-crystalline. They include several widely distributed polymers (such as PE, PA, PET and PP [5,6]), as well as high technology polymers (aromatic PEEK or PEKK for instance [7]). For these materials, crystallization, which occurs during the forming process, and especially during the solidification phase, plays a critical role in the final quality of the part. For instance, in the case of injection moulding, crystallization determines the liquid/solid transition and thus rules the ejection criteria. In addition, crystallization induces shrinkage and residual stresses.
Therefore, polymer-forming simulation softwares must accurately predict crystallization effects, especially during the cooling phase. An accurate quantitative prediction is a prerequisite for predicting the solidification (or gel) time.
The crystallization of polymer macromolecules depends on both temperature and flow history. A macroscopic homogenized approach consists of defining a continuous crystallinity variable x, which evolves from 0 for an amorphous material to a maximum crystallinity x max , specific to the semi-crystalline polymer considered. The degree of crystallization α is also defined as the degree of advancement of the transformation and thus varies from 0 (when x=0) to 1 (when x=x max ). The evolution of α is ruled by the crystallization kinetics law [8]. This paper focuses on quiescent thermal crystallization. In fact, temperature appears to be the prevailing driver of crystallization in several forming processes, in which flow-induced crystallization can be neglected.

Avrami Integral Form
The most common thermal crystallization models are derived from the Avrami model [9]. Under Avrami assumptions (uniform germ location, instantaneous growth after activation, given crystal growth morphology [10]), isothermal crystallization evolution is written: where n is the Avrami index and K av the Avrami crystallization kinetics function. In the seventies, Nakamura et al. [11] extended the Avrami model to account for non-isothermal crystallization. In its initial form, the degree of crystallization was given as the integral form: where K ( ) T is the Nakamura kinetics crystallization function which is related to the Avrami isothermal kinetics . This expression was also obtained by Billon et al. [12,8]. The Nakamura model proved efficient for modelling the crystallization of different polymer systems [13,14,15].

Differential Form
Patel & Spruiell [5] suggested using the differential form of equation (1) which is written The Nakamura function G is plotted over [0,1] in Figure 1. This differential form, a first-order non-linear ordinary differential equation, is more natural to implement and solve numerically, especially in a multiphysical framework [16,17,18,19,20,21,22,23]. The initial state of the material (amorphous or semi-crystalline) should be given as an initial condition for α.
Even though several polymers exhibit more complex crystallization behaviour, the Nakamura law (2) very often appears in the kinetics. This is the case, for instance, with kinetics accounting for secondary crystallization, in which two coupled Nakamura laws can be used [24,25].

Artificial Germination
In its differential form, because G ( ) 0 =0, an initial fully amorphous state α ( ) t=0 =0 would never induce crystallization. Thus, it is common to consider an artificial initial germination where ε is a small value. To the author's knowledge, its influence has never been fully investigated in the literature and will be discussed below.

Time Integration
Equation 2 is highly non-linear. In the case in which the temperature T ( ) t (and thus K ( ) there is no analytical solution to this equation. Therefore, a numerical integration is required.
The time and degree of crystallization are discretized. In this short communication, a standard first-order incremental integration scheme is considered. At an instant t n , the degree of crystallization α n is computed from the previous value α n−1 .
Two extreme and very general integration schemes are considered: 1. An explicit forward Euler integration scheme which gives rise to a constraint on the length of the time step to ensure stability. α n is obtained as δt being the time step.
2. An implicit backward Euler formulation, which is stable but requires a non-linear resolution at each time step. α n is obtained by zeroing the residual Equation (4) is usually zeroed with an iterative method, which may require the computation of the derivative of dG/dα. Difficulty. During iterations, the Nakamura function G and its derivative could be evaluated for various values of α. Even non-physical values ( α∉ [ ] 0,1 ) could be tested. In these cases the residual evaluation should return a value that allows further convergence.  (2)) is singular when α=0 and when α=1 (vertical tangent, represented by the red arrows). In the proposed modified extrapolation, a thresholding is performed below α min , and a linear extrapolation above α max .

Modified Nakamura Function
The Nakamura function G is modified in order to ensure robust numerical evaluation. As depicted in Figure 1, the original Nakamura function G exhibits two singular points with infinite derivatives for α=0 and α=1.
These infinite derivatives may prevent convergence of the nonlinear solving of equation 4, especially with gradient methods. Two extrapolations are proposed on the function G.
For values of α below a minimum threshold α min , G ( ) α is truncated to G(α min ). This ensures: (i) an artificial germination for the start of the initial crystallization,

Results and Discussion
In the following, in order to focus on the efficiency and robustness of numerical integration schemes, the study is performed with a constant The reader should keep in mind that in real industrial cases, the temperature field (and therefore K ( ) T ) depends on space and time. Thus, equation (2) should be integrated in a global multiphysical integration scheme. For such applications, the reader may refer to references [16,17,19,20,21,22].
Considering this degenerated case, with K ( ) T =1 , the Nakamura equation (2) has an explicit analytical solution given by equation (1). This is helpful to study the convergence of the numerical schemes: in particular, the crystallization half-time for which α reaches 0.5 is known:

Lower Truncation α min
The effect of the lower truncation α min , which accounts for the artificial germination, is first investigated. Two classic values of the Avrami index n=2 and n=3 are considered.  If α min is too high, the artificial germination is exaggerated and the crystallization time is underestimated. This parametric study suggests that the threshold value α min should be below 10 −6 for Avrami index n=2, and below 10 −8 for Avrami index n=3.

Upper Truncation α max
As discussed in section 2.2, the time integration scheme may result in evaluating the Nakamura function G for values of α above 1. This may occur in the cases of explicit formulations with large time steps or implicit formulations during the non-linear resolution (4). The upper linearization of G ensures robustness in these cases.
With an explicit time integration with time steps dt=0.2 s, Figure 3 shows that a slight overshoot occurs. Nevertheless, the extrapolation of the Nakamura function above α max forces α to approach 1 in further increments. Thus, the proposed extrapolation enables a robust time integration. Figure 4 shows the convergence of the crystallization half-time t 1 2 versus the size of the time step δt considered. It is advisable to use time steps no larger than [5.10 −3 ]s to ensure accurate kinetics prediction 1 . In an industrial situation, faster time integration schemes (such as Crank Nicholson or Runge-Kutta) might also be used. Figure 5 is the convergence plot representing the relative error for both schemes versus inverse of time step. The classical first order convergence is 1 In the general case in which the Nakamura function K ( ) T differs from unity, this time step recommendation should be divided by the function value. obtained with explicit or implicit scheme. Classical results in terms of stability can also be recovered. For instance, the resolution becomes instable for large time steps (dt>0.3s) with the explicit scheme.

Conclusion
The crystallization of polymer macro-molecules is usually modelled by defining a macroscopic degree of crystallization α, which varies between 0 and 1. Classic modeling of crystallization kinetics is based on the Nakamura equation. This paper investigated the accuracy and robustness of the numerical time integration of Nakamura crystallization kinetics. The main contributions of this short communication are the following: For robust integration schemes, the Nakamura function should be truncated (i) for small and negative values of α and (ii) for values of α in the vicinity and above 1 Guidelines concerning these truncation values are given. They are fairly restrictive (α min <10 −6 ).
Guidelines for time stepping, in the case of standard forward and backward Euler schemes, are also provided.