Simulation framework for crystallization in melt flows of semi-crystalline polymers based on phenomenological models

Polymer components are shaped mostly out of the molten state. As in the case of semi-crystalline polymers, crystallization can be suppressed by shock cooling, thermal process design allows to influence the solid bodies properties. A simulation approach that enables to predict these properties based on a forecast of crystallinity is presented in this paper. The main effects to consider and possibilities of modeling and simulation are discussed. A detailed description of how to create an experimental foundation using dynamic scanning calorimetry (DSC) and a rheometer is provided. Suppression of crystallization is modeled by a novel phenomenological approach, based on data over a large band of cooling rates. Special focus is put on parameter identification and extension of insufficient DSC data. The mechanical behavior is modeled using a weighted approach based on a nonlinear-thermoviscoelastic model for the molten state and a highly viscous Newtonian model for the solid state. Parameterization of both models is highlighted. An implementation in OpenFOAM is documented, emphasizing specific methods that were applied. Results of simulations for a simplified profile extrusion and injection molding case are presented. Basic relationships are forecasted correctly by the method, and important findings are presented for both processes.


Introduction and classification of the method
Typical manufacturing processes for polymer components are based on three basic steps: melting, conditioning and shaping. For instance, a granular might be melted in the first few elements of an extruder, after which it is degassed and fed into a profile extrusion tool thereafter, see [1] and [2]. Solidification is initiated in the final step to form a solid body. By thermal control of the cooling process, opportunities to influence the properties of the solid are given. This is most relevant for semi-crystalline polymers, given that their crystallinity is highly dependent on the applied rate of cooling. Crystallinity, e.g., influences density [3], opacity [4] and mechanical properties [5]. Therefore, a simulation method that allows predicting crystallinity is of great interest.
The underlying mechanism is that, during non-isothermal crystallization, the applied cooling rate influences crystal growth. In the case of low cooling rates, large crystals develop, leading to high crystallinity. For high cooling rates, many small crystals can be found, resulting in a large count of impingements that reduce crystallinity, see, e.g., [6]. Furthermore, it is possible to completely suppress crystallization for even higher cooling rates. The standard procedure for investigating these effects is to perform DSC [7] scans. Only recently, commercial equipment that allows to investigate the latter effect at such high cooling rates was made available, see [8] and [9]. This makes possible a new way to formulate phenomenological crystallization models. In combination with modeling the influence of crystallization on a flow process, this represents the core of this paper.

Crystallization modeling
For modeling polymer crystallization, three common approaches exist. Firstly, which is probably the most popular, the theory of volume-averaged expanding spheres. It dates back to the works of Kolmogorov, Avrami and Evans [10][11][12] summarized, e.g., in the book of Janeschitz-Kriegl [13]. The most used formulation is based on four differential equations, known as Schneiders rate equations [14]. Given that it enables the inclusion of processes starting from the scale of nucleation, it is possible to consider a wide variety of effects that influence crystallization. However, it is not intuitive to use, since a direct connection to observations from standard DSC scans is difficult. Secondly, there is an approach that only uses one differential equation. It is known as the model of Nakamura [15] and Ziabicki [16], which can be directly fitted to standard DSC experiments. Therefore, it might be regarded as more intuitive. Finally, there are phase field methods [17]. A well-known field of their application is crystallization of metals [18]. They are applied rarely for polymers, see, e.g., [19]. Equations are usually derived from the Clausius-Duhem inequality and phenomena are described on a micro-scale. This allows simulating the growth of actual crystals. However, usage in problems with technical scale is not possible without applying a homogenization technique. For all models that are named above, the implementation of findings from DSC experiments at a fixed cooling rate is possible. An adaptation to experiments over a wide range of cooling rates, which are investigable with modern equipment, is not known. For this reason, the present work's approach is derived directly from such experiments.

Rheological modeling
Even though crystallization modeling is important, many more effects need to be considered. Considering that polymer melts are complex fluids, standard approaches, e.g., the usage of a generalized Newtonian model [20], is not sufficient. Common effects to be modeled are shear thinning, normal stress differences, strain hardening and stress relaxation, which are additionally temperature-dependent, see [21] and [22]. For this purpose, the usage of Maxwell type of equations has been established. Common models applied are the Giesekus [23], Phan-Thien and Tanner [24,25] or finitely extensible nonlinear elastic (FENE) [26] models. So as to make constitutive equations temperature-dependent, time-temperature superposition is used widely. Its concept is described in [27], the application to constitutive equations in [28]. Furthermore, polymer melts are known to be highly dissipative [29], which requires consideration in the energy equation. In the case of Maxwell type of equations, a special treatment is developed in [28]. When considering crystallization as a fluid-solid transition, it becomes clear that it is a flow dominating effect. In order to model its impact on mechanical behavior, it is common practice to use formulations causing a strong increase in viscosity, summarized in [30].

Classification of the present work
As described in 1.1 and 1.2, simulating crystallization in polymer melt flows touches many disciplines, and in each, a large variety of models are available. Combining these to create a framework for predicting crystallinity is a challenging task, especially because an experimental foundation for each approach has to be available, since material behavior is very sensitive to constituents and molecular structure of the polymer. Presumably, hence, not many authors presented a framework suitable for this purpose. In 2014, Zhao [31] presented a solver based on Schneiders rate equations. The melt model of this work is based on a FENE approach, the solid regions are based on an orientation-based stress tensor. An influence of deformation rates on crystallization is considered, and the parameters used originate exclusively from the literature. The framework is based on the Finite Element Method and simulations are performed for a profile extrusion tool. In 2016, Zhao and Mu [32] presented a second paper in which they changed the solver to a Phan-Thien Tanner model for the melt and a viscous formulation for the solid. In the same year, Spina [33] published a paper on simulating injection molding using a two-phase solver. For crystallization, a two-scale model adding structure simulations to the model of Nakamura [15] is used. Both fluid and solid states are modeled as non-isothermal purely viscous fluids. An influence of deformation rates on crystallization is not considered, and the models are chosen in reference to experiments. A finite element-based solver was used as well.
Approaches applied in the above-mentioned works are used in the present work as well. The greatest difference is the crystallization model, as it is derived directly from experiments that investigate the suppression of crystallization. Therefore, crystallization enthalpy is not a constant, and it is possible to investigate the suppression of crystallization numerically. In terms of fluid modeling, viscoelastic constitutive equations are utilized, but the description of the solid is based on a limit value of viscosity. For testing the method, simulations are performed for a profile extrusion and injection molding case. In contrast to the works named above, a finite volume method-based solver is used.
The benefit of the presented crystallization model is that it can be implemented in any CFD code straightforwardly. It is implemented in the open-source software OpenFOAM in the present work. For demonstration, because of its great relevance in research and technology, the framework is adapted to isotactic polypropylene. The specific material is Sabic PP 575P in granular form, an unfilled multipurpose polymer for injection molding.

Experimental foundation
The simulation framework of this paper uses models to describe the crystallization and rheology of iPP. In this chapter, the general methodology to obtain an experimental foundation is described. Particular focus is placed on general dependencies and common problems that occur when using standard equipment.

Thermal analysis
In order to gain information about the crystallization behavior, investigations using differential scanning calorimetry (DSC, [7,34]) were performed. The equipment used are the TA-Instruments DSC Q1000 and PerkinElmer Pyrus DSC 6, which allow performing measurements at moderate cooling rates. Therefore, it is not possible to investigate the suppression of crystallization using these machines. In short, a DSC is equipment able to detect any heat flux emitted or absorbed by a sample. Experiments are performed at constant heating or cooling rateθ over a certain temperature range. In case a sample shows constant heat capacity c, the resulting DSC signal can be described asq DSC (θ ) = −cθ +q(θ,θ). (1) So as to give an insight into the melting and crystallization behavior of iPP, measurements were performed between 0 and 250 • C at three temperature rates. The results are presented in Fig. 1 by plotting the apparent heat capacity c ap =q DSC /θ. Arrows annotate the experiment's thermal direction.
Considering the peaks in this plot, it shows that, e.g., for 2.5 K min −1 , the sample melts at around 165 • C and crystallizes at 130 • C. From this, an approximate discrepancy of 35 K results, which is detectable for all cooling rates presented in Fig. 1. This finding means, just crystallized regions will not melt because of a small rise in temperature. Regarding crystallization processes in technical devices that are designed particularly for cooling, this means crystallization can be considered to be an irreversible process, as such devices do not allow the temperature to increase. It is a valuable effect for modeling considering that this demands a model formulation allowing only progress of crystallization. In order to investigate the crystallization behavior further, a greater range of cooling rates was investigated. However, since only the latent heat signalq(θ,θ) is of interest, the caloric signal is subtracted based on a polynomial fit. This leads to graphs shown in Fig. 2, already allowing to determine general dependencies for crystallization.
Firstly, the latent heat release causes a Gaussian bell-shaped heat flow signal. Secondly, with an increasing cooling rate, peak magnitude and width increase, additionally the peak temperature is shifted to lower temperatures. Furthermore, with respect to the complete suppression of crystallization, which means that no latent heat is released, no such tendencies can be observed in the presented data. Consequently, the cooling rate has to be increased further. Results of such experiments are presented in Fig. 3.
It is shown that the general relations already discovered regarding Fig. 2 can also be discovered for higher cooling rates. However, in case of the highest rates presented, initially visible for 25 K min −1 , the symmetry of the bells is lost. This effect is not to be mistaken as material behavior. Rather, it is caused by the sample  mass being too large to sustain uniform crystallization throughout the sample. A well-known effect to DSC users that is solvable by choosing a lower sample mass. However, since the sample mass that was used can already be considered small in the context of standard DSC equipment, an impasse is reached. Depending on the DSCs specifications, any user will face this problem at some point. The solution is to extend measurements with data obtained in nanocalorimetry. The latter is a new field in thermal analysis which received increased attention as novel commercial equipment is available, e.g., presented in [8]. With this equipment, even for cooling rates at which there is no sign of crystallization anymore, data can be gathered in great detail, as the named publication shows. Such a machine could not be used in the present work. For this reason, an approach that allows extending experimental data was developed, based on a special method of evaluation that relates to the modeling approach chosen in the present paper. Therefore, the topic of crystallization is reconsidered in Chapter 3.1 and the current chapter continues with the influence of crystallization on rheology.

Rheological analysis
Considering that the developed simulation approach also covers the influence of crystallization on rheology, cooling experiments were performed in a rotational rheometer, e.g., described in [35]. The rheometer used is a TA-Instruments AR-G2 equipped with a testing chamber that makes possible investigations with a parallelplate system in a nitrogen atmosphere. The plates diameter is 25 mm and the measuring gap was set to 600 µm which corresponds to an approximate sample mass 280 mg. At a constant shear rate ofγ = 0.01 s −1 , cooling was initiated for different rates starting from 200 • C, for which the results are shown in Fig. 4. DSC curves plotted for comparison were recorded using the Q1000 with a 7.74 mg sample.
The investigations show that, as soon as crystallization sets in, a singularity occurs in viscosity. However, compared to a DSC, the sample mass is much higher here. Consequently, uniform crystallization throughout the sample cannot be expected. With regard to the setup that was used in the present work, it was even visible to the naked eye. Crystallization started at the outer surface, moving inward over time. Consequently, those results should not be used to connect the state of crystallization with viscosity. The only effect that can be drawn from these experiments is that, already at the onset of crystallization, a great increase of viscosity can be expected.

Modeling
Key of this work's crystallization model is the fluid-solid indicator χ , which ranges from χ = 0 for fluid to χ = 1 for solid. Based on its evolution, which is described by the crystallization model, terms in the balance equations are changed. In the energy equation two additional source terms are present. The first is S fluid , which considers dissipation. It is blended out during crystallization; hence, only heat conduction is considered in the solid phase. The second source term considers the latent heat release during crystallization. It uses the crystallization enthalpy h cryst as a function of the material time derivative of temperatureθ = Dθ Dt . Stress is modeled by blending fluid stress T fluid and solid stress T solid linearly based on the phase indicator. This approach is considered in the momentum equation as follows: Changes of density ρ, heat capacity c and heat conductivity k are not considered in this work, since continuous manufacturing processes are prioritized. The used parameters are documented in Table 1. In case of discontinuous processes like a complete injection molding cycle which includes the holding stage should be modeled, dependencies given in [36,37] need to be included.

Crystallization
In the present work, crystallization is described by a single differential equation, as is the case in the model of Nakamura [15] and Ziabicki [16]. However, instead of deriving growth rates from experiments, which are then used to describe crystal growth volume averaged, a more phenomenological approach is derived. The latter is based on the evolution of the fluid-solid indicator, which can be directly gathered from DSC experiments.

Phase transition
So as to describe the phase transition, a function inspired by formulations of Ziabaicki in [38] is used, namely the normalized temperature integral of the crystallization rate, used in his model to drive crystallization, rewritten as To make it cooling rate dependent, parameter functions are considered. For identifying these, mid-temperature θ m and half-width of the crystallization rage D 1/2 had to be determined for a great range of cooling rates. For this purpose, the latent heat flux of each measurement needs to be integrated over temperature and taken in reference to the applied cooling rateθ exp . The result corresponds to the enthalpy release, highlighted in the following equation: The way in which experimental data are transformed by the application of this integration can be seen in Fig. 5. It was carried out numerically using the trapezoidal rule. The according value of crystallization enthalpy can be determined easily, since it is the value reached after the heat flow vanishes. This allows to experimentally determine h cryst (θ), used in Eq. (2). For finding the remaining parameter functions, it is helpful to know that Therefore, after h cryst is determined for the applied cooling rate, the mid-temperature can be determined as θ m = θ( h cryst /2). The half-width D 1/2 is then identified using least squares. For each cooling rate, this results in one data point for each parameter, as shown in Fig. 6. All experiments that were performed for the present work, including those presented in Chapter 2.1, are plotted in red. Comparing those to data points for h cryst that were determined for measurements presented in literature shows that, indeed, the self-investigated cooling rate range is far away from the range at which crystallization is suppressed. Therefore, as already described at the end of Chapter 2.1, the aim is to extrapolate data. The presented evaluation strategy is suited for this purpose excellently, as Fig. 6 shows.
In material science, h cryst (θ) and θ m (θ) are well investigated for iPP. However, it is not common practice to present bell curves for a large range of cooling rates. Therefore, the choice of publications that enable identifying D 1/2 (θ) is smaller. Enough sources could be found nevertheless, since iPP is a material of great interest, see the bottom plot of Fig. 6.
Prior to introducing the crystallization model, key points that can be gathered from the data must be identified. The data for h cryst (θ) shows that, at around 250 K s −1 , suppression of crystallization sets on to a great extent. Additionally, with an increasing cooling rate, the mid-temperature θ m decreases steadily. Finally, D 1/2 is nearly constant below 0.1 K s −1 and starts to increase for higher rates. On the basis of data determined for [8] and [39], it can be suggested, that there is a limit value for high cooling rates for this parameter. However, the investigated cooling rate range was not sufficient for a certain statement.
Parameter functions were formulated based on these findings. The according parameters are documented in Table 2. Some of these were identified manually, mentioned in reference to each approach. All others were identified using least squares. So as to model the crystallization enthalpy, an approach by Schawe [8] was used. It is based on a function for generic crystallinity g(θ) and retardation v(θ). However, in the way g(θ) is proposed in [8], it will turn negative for vanishing cooling rates. For this reason, it was modified in the present work. The resulting model for crystallization enthalpy is Interpretable parameters are h 0 = h cryst (θ = 0), which should be read from a graph using linear scale iṅ θ . This value represents the maximum possible crystallization enthalpy release, belonging to the highest value of crystallinity for the material. The value of v 2 is the cooling rate, and strong suppression of crystallization sets in.
To describe the mid-temperature of crystallization, a power law in the form of is chosen. Just as h 0 , the zero value θ m,0 = θ m (θ = 0) can be determined in a linear plot of the experimental data. For D 1/2 , since the data from Fig. 6 allow the assumption of a limit value, a Cross [42] function is used. D 1/2,0 is the value of D 1/2 at low cooling rates, D 1/2,∞ is the limit value. The reciprocal value of the cooling rate, at which D 1/2 starts to increase, is the approximate value of k D . A critical property of those functions is that they are not usable for positive temperature rates. In the simulations performed for the present work, this state will not occur, as sufficient cooling is provided. Nevertheless, the limiterθ = min Dθ Dt , −1e-6 (10) was used for implementation in the code.

Evolution equation and evaluation of crystallinity
In order to consider the phase transition in simulations, an evolution equation for χ needs to be formulated. It is derived calculating the total differential In contrast to a phase-field model, which would describe the progress of the fluid-solid interface, this equation describes the progress of the solid's volume fraction. The resulting interface is homogenized and diffuse, hence the application is limited to quasi-static crystallization processes. Considering that the fluid-solid indicator is bounded between χ = 0 and χ = 1, for evaluation a discrete interface can be constructed at χ = 0.5. However, it is not interpretable as a sharp, physical fluid-solid interface. In Eq. (11), the term includingθ is neglected. This is possible for processes that include strongly cooled boundaries which do not allow the latent heat to dominate and, e.g., cause local heating. In the remaining term, the temperature derivative is calculated analytically, whereas the material time derivativeθ is calculated numerically in the simulation. However, the formulation of the equation is not complete yet, as one important property is missing. In reference to Fig. 1, it was mentioned that the crystallization model should not allow melting. In terms of the evolution equation, this means only positive right sides are allowed. This property is implemented by a limiting function, leading to As described in reference to Eq. (2), the rate of change for χ is used in the latent heat source. Therefore, it is possible to derive the heat released by integration, denoted Q. Solving the additional transport equation (13) adds this feature to the solver. It allows to derive the relative crystallinity c r by taking the final value of Q in reference to the maximum possible value of crystallization enthalpy h 0 from Eq. (7). Local evaluation is enabled by The value obtained is relative to the material's specific maximum possible crystallinity. For iPP it can be calculated by taking h 0 in reference to the crystallization enthalpy of pure crystals. According to [43], this value is 207 kJ kg −1 . Hence, the value of c r = 1 would represent an absolute crystallinity of 56 %.

Melt rheology
Polymer melts are considered as viscoelastic fluids. The method to describe them utilizes the elastic viscous split stress (EVSS, [44]) approach. Accordingly, extra stresses are modeled as the sum of polymer (p) and solvent (s) stresses, and hence, Solvent stresses are proportional to a residual Newtonian solvent viscosity. Consequently, they modeled as based on D = 1 2 (L + L T ), the symmetric part of the velocity gradient L = ∇v. Polymer stresses are described based on a relaxation time λ, a non-linear extension term Q(T p )T p and the polymer viscosity η p . They are considered in a Maxwell type of transport equation in which the upper convected Oldroyd derivative occurs. Based on these definitions, it is possible to describe the solvent contribution by which is determinable from steady-state shear experiments. Those are measurements performed in rheometers, that cover a certain range of the shear rateγ . It is common practice to determine zero-and infinite-viscosity from them by for evaluation. Since in EVSS, η ∞ = η s and η 0 = η p + η s , the value of β can be calculated with no additional effort. However, it is a challenging task for most polymers to determine η s as the shear rates needed are too large for standard equipment [45]. Measurements presented in [46] and [47] show that for iPP, β = 0 can be chosen, as the values measured for η s are 4 magnitudes below η 0 . Consequently, the solvent viscosity is set to η s = 0. The choice of an appropriate nonlinear extension depends on the elongation behavior of the melt. In [48][49][50][51], it is shown that for moderate strain rates, iPP behaves weakly strain hardening, with a strain hardening ratio of 1.2. Furthermore, for rates above it behaves strain softening. These effects are implemented using the exponential extension of Phan-Thien and Tanner (EPTT, [25]) with α = 0.5. So as to make the transport equation temperaturedependent, time-temperature superposition in the sense of [28] is applied. For this purpose, the Arrhenius equation is used and a temperature dependency of polymer viscosity and relaxation time is introduced by Its parameters were identified by mastering results of steady-state shear experiments performed at different temperatures, shown in Fig. 7. Given that it is generally not possible to fit the EPTT model over the entire shear rate band, it was fitted to the Newtonian range and the range of shear thinning. Additionally, the cooling experiments performed for investigating crystallization, already presented in Fig. 4, were used. The identified parameters, used for the master curve that is presented Fig. 8, are documented in Table 3.
To consider dissipation in viscoelastic fluids, standard treatment is not sufficient, as is described in [28]. An additional source term ought to be considered, using a weighted approach for the dissipation source term of Eq. (2). Only the implementation of the first term is known to produce negative dissipation; hence, the second term is weighted stronger by using κ = 0.25.

Solidification
Modeling the impact of crystallization on mechanical behavior is generally conducted by manipulation of viscosity, best summarized in [30]. The present work implements this behavior by the usage of a finite viscosity value η solid . Stresses in the solid phase are then described by the Newtonian model In order to choose a reasonable value, a 1D pressure-driven channel flow was solved analytically for a spatially changing viscosity. A similar problem is, e.g., treated in [52], but for a different purpose. As can be seen in Fig. 9, in the top half of the channel, viscosity is set to η 1 , in the bottom half it is set to η 2 . The two fluids are therefore Newtonian; at the interface the coupling condition prescribes that velocities and stresses are equal. A variation of parameters shows how increasing viscosity in the lower part decreases the interface velocity, c.f. Fig. 9. The following was discovered: in reference to the maximum velocity u max = max[u(y)], a ratio

Simulation approach
The solver in which all models are implemented is composed in foam-extend-4.0, a fork of the open-source CFD library OpenFOAM. For pressure correction, it uses the block coupled method described in [53]. So as to enable a stable treatment of viscoelastic models, the log-conformation reformulation (LCR) [54] is implemented based on the methods described in [55] and [56]. The convective and diffusive terms are treated implicitly in all transport equations. Euler time integration is applied, and time step size is controlled based on the Courant number, which is kept below 0.1. Convection terms are discretized using the upwind method, central differences are applied for diffusion terms. The top-level algorithmic structure is documented in Fig. 10.
Updating the Arrhenius functions means that θ is inserted in Eq. (21). For full pressure correction, a block coupled pressure correction loop, not displayed in Fig. 10, is repeated until the residuals for p and v fall below 1e−5. In the stress correction loop, which updates T p , Eq. (20) is inserted in Eq. (17), of which the latter is solved afterward. Only if the initial residuals of pressure correction and stress update are below 1e−5, the algorithm proceeds. This corresponds to a strong pressure-velocity-stress coupling. The next step solves Eq. (12) for the progress of crystallization, which is used afterward to integrate the latent heat by solving Eq. (13). Finally, the energy equation is solved, since all its source terms depended on fields computed beforehand, c.f. Eq. (2).
Since computational effort is always of importance, the expense of each step is discussed in the final part of this chapter for a 2D case. It is best explained with the number of unknowns (N ), which, for the finite volume method, corresponds to the number of cells. As common for incompressible flows solvers, the most expensive operation is to obtain simultaneous fulfillment of the momentum and continuity equation by pressure  In order to validate and verify this solver thoroughly, a three-staged approach was used, of which the results are published in [57]. In stage one, only the flow section of the code is validated using a lid-driven cavity and startup channel flow. Stage two covers the crystallization section of the code. It is validated with simulations for a virtual DSC at different cooling rates and verified on the Stefan problem. In stage three, the overall code with both sections activated is verified on crystallizing 1D and 2D flows. In the reference named, also a mesh study was carried out for a crystallizing lid-driven cavity flow.

Examples of application
The use case of the method developed is intended to predict crystallinity. In order to provide two examples of application, simulations are performed for a flow channel inspired by a real tool presented in [58]. As shown in Fig. 11, the process is divided in two sections. Table 4 Boundary conditions used in the test channel. Note that the cooling section's wall boundary conditions are named afterward. For profile extrusion in Table 5, for injection molding in Table 6 Boundary Velocity Pressure Temperature The front section exists to obtain uniform flow. Crystallization is induced by cooling in the back section. The dominating dimensionless number for this problem is the Peclet number formulated using channel height h and inlet velocity U . In a tool with these dimensions, process velocities in the low centimeter per second range are realistic. Therefore, when considering 1 cm s −1 as a characteristic value Pe ∼ 520, meaning that the process is convection-dominated. The domain is meshed using 50 uniform cells over the height and 700 over the length. By choosing certain velocity boundary conditions in the cooling section, an exemplary setup representing profile extrusion or injection molding can be set. A special temperature boundary condition is used in both cases. It prevents having to use constant wall temperature, which would result in a jump of temperatures between the walls of conditioning and cooling section. The motivation is to mimic water cooling based on a heat transfer coefficient α w = k wall (k d wall ) , which is derived from the law of static thermal transmittance. Accordingly, the heat flux exiting the polymer flow is equal to the heat flow entering the water, described by Heat conductivity of the wall material k wall , wall thickness d wall and coolant temperature θ c are considered in this equation. By rearranging Eq. (26), the inhomogeneous Neumann boundary condition ( ∂θ ∂ y ) is derived. When assuming a steel wall (k wall = 50 W m −1 K −1 ) with a thickness of d wall = 2.5 cm, α w ≈ 10.000 m −1 . The coolant temperature is set to θ c = 10 • C. All remaining boundary conditions are specified in Table 4.
For velocity, a block profile is set at the inlet to ensure a constant volume flux with minimum effort. Development of a parabola-shaped velocity profile will take place in the first few millimeters of the conditioning section. At its walls, a no-slip boundary condition is set. Since the melt is still fluid in this range, no wall slip is expected here. At the cooling sections outlet, the velocity boundary condition for an undisturbed outlet is used. For pressure, the reference value is set at the outlet. Since the flow is incompressible and therefore the total pressure value is not relevant, it was set to zero. Given that pressure drops in flow direction, this ensures only positive values. At all other boundaries, pressure will develop during the simulation. The inlet temperature is set to 200 • C because it is close to the crystallization range, see Fig. 1. This way, crystallization will take place in a confined space after the beginning of the cooling section. For fields not named in this table, homogeneous Neumann boundary conditions are set. This means boundary values are computed regarding the corresponding cell values. As initial condition, the fully molten, static and stress-free state at 200 • C is assumed.

Profile extrusion
So as to simulate a profile extrusion case, the boundary conditions named in Table 5 are set in the cooling section. The top wall is water-cooled, and the bottom wall is adiabatic. Such a setting would be chosen for, e.g., hollow profile extrusion, see [58]. It is well known that, due to crystallization in extrusion tools, wall slip occurs. Therefore, the velocity boundary condition is linked to the phase indicator. It is used to interpolate between no-slip and full-slip linearly, which can be realized by using the cell centroids horizontal velocity u P of the boundary cell.
The steady-state result for U = 0.1 mms −1 is depicted in Fig. 12. In this case, heat diffusion is strong enough to let the melt fully crystallize inside the computational domain. For the steady-state, stream lines represent path lines. This allows to determine that, with regard to particles closer to the top wall, the traveled path during crystallization is shorter than for those in the bottom wall region. Furthermore, at the cooling section's inlet, the melt stream is slightly deflected downward by the premature crystallization at the top wall, best visible for the horizontal velocity u. Due to crystallization, slip is initiated on both walls. Therefore, the solid leaves as a block with uniform velocity. The temperature field is typical for a flow that is cooled on one side. Due to the low extrusion velocity, penetration depth of the cooling is high. Increasing the extrusion velocity stepwise allows the discovery of two relationships, see Fig. 13. Firstly, increasing U leads to the melt not crystallizing fully inside the domain. The thickness of the crystallized layer is decreased strongly, as the penetration depth of the cooling is decreased. Secondly, due to an increased process velocity, the cooling rates are increased, visible in the lower relative crystallinity. It is caused by convective cooling, which is part of the material time derivative of temperature.
This opens up the possibility for interesting studies. For example, at the upper wall, relative crystallinity does not develop over length. Full crystallization is reached at the cooling section's inlet, from which can be concluded that the outer layer's crystallinity is only influenced by the thermal transition from conditioning to cooling section. Therefore, if, e.g., higher crystallinity is important, a smoother transition has to be realized.

Injection molding
The cooling section's boundary conditions for modeling an injection molding case are named in Table 6. Just as in a typical mold, cooling occurs from all directions, in consequence both walls are cooled. Considering that there are smaller pressure gradients in injection molding, usage of no-slip boundary conditions is sufficient.
In light of the fact that extraction of heat by two boundaries is more effective, a higher inlet velocity is chosen. Focus of study is the interaction of molten and solidified regions. In Fig. 14, the steady-state result for As a result, velocity is increased in the center. The temperature distribution shows that the core remains almost at the inlet value. Dissipation contributes to thermal isolation of the core. It is proportional to stresses and velocity gradients, c.f. Eq. (23) and therefore has its maximum in the regions close to the interface. Shear and normal stress develops as expected for a contraction flow.
Analogous to the profile extrusion case, a strong influence of process velocity on the shape of solidified regions can be observed, see Fig. 15. However, crystallinity in proximity of the mold walls is not influenced noticeably. This is a result of the melt not slipping at the walls and a parabolic velocity profile occurring for this range of process velocities, c.f. Fig. 14. Hence, cooling is diffusion-dominated in the investigated cases. Regarding this observation, the distribution of crystallinity over the thickness of solidified regions is plausible. Just as for a static cooling case, at process beginning, the temperature gradient at the wall is high, leading to the largest cooling rates in the process. This results in the lowest crystallinity. In progress of cooling the temperature gradient smoothens out, leading to smaller cooling rates, causing an increased crystallinity. Therefore, the largest crystallinity is always located at the fluid-solid interface.
The conclusion of this study shows that, if wall crystallinity should be decreased, much larger injection velocities are needed. This leads to block-shaped velocity profiles, as they are common for shear-thinning fluids, which cause a larger temperature gradient at the walls. To find the right injection velocity for a desired crystallinity, the method presented in this paper can be used.

Summary and outlook
This work presents a novel approach to consider cooling rate dependent crystallization of semi-crystalline polymers in CFD. It is demonstrated by its application to isotactic polypropylene. Opportunities to gather data for crystallization modeling from thermal analysis and rheometric investigations were presented. In order to investigate crystallization, DSC measurements need to be performed over a wide range of cooling rates. However, with standard equipment, high cooling rates cannot be reached. Investigations performed with a rotational rheometer showed that the sample mass in such a measuring system is too high to find a relationship between the state of solidification and a rise in viscosity. Nevertheless, it is undisputed that, already at the onset of crystallization, the viscosity increases strongly. The crystallization model developed in this work is an evolution equation, set up for a fluid-solid indicator, motivated by a Gaussian error function that can be fitted to the DSC data. These data are gathered from own experiments, complemented by data from publications that  15 Influence of the injection velocity on the final relative crystallinity use nanocalorimetry. Stress in the molten state is modeled using the exponential Phan-Thien Tanner model made temperature dependent by using an Arrhenius function. The solids' stress is described with a Newtonian model using a solid viscosity 1000 times larger than a reference viscosity for the molten state. All models were implemented in foam-extend-4.0 into a solver that was composed to use block-coupled pressure correction and the log-conformation reformulation. So as to demonstrate the capability of the approach, results of simulations for a simplified profile extrusion and injection molding case are presented. They showed that basic relationships are forecasted correctly by the method developed. It was discovered that surface crystallinity of profile extruded products can be influenced by the thermal transitions at the walls of a shaping tool. Furthermore, it is possible to influence surface crystallinity with the injection velocity in injection molding.
Further development is possible after validating the method with experimental data. This could be the crystallinity distribution in an extruded profile. Gaining this information is very challenging, since it requires investigating small samples from the specimen optically or calorimetrically. Those need to be taken at ample different locations and the process of cutting them out will probably change the crystal structure. However, this methodology will provide a first insight that is valuable for the progress of this work. There are still many ways to adapt the method to experimental findings. In development of the crystallization model, e.g., mechanical load was neglected. There is, however, an influence of shear, elongation and pressure which is documented extensively in [13]. Model extension based on an experimental foundation is possible. One approach could be to extend Eq. (12) by a dependency on shear and elongation rates. Therefore, by using e.g., a difference in the impact of shear and elongation can be considered. It is a well-known effect with great relevance in manufacturing, see [59]. Preliminary work considering shear has been documented in [57].
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Funding Open Access funding enabled and organized by Projekt DEAL.