Modeling Injection Molding of High-Density Polyethylene with Crystallization in Open-Source Software

This work investigates crystallization modeling by modifying an open-source computational fluid dynamics code OpenFOAM. The crystallization behavior of high-density polyethylene (HDPE) is implemented according to theoretical and experimental literature. A number of physical interdependencies are included. The cavity is modeled as deformable. The heat transfer coefficient in the thermal contact towards the mold depends on contact pressure. The thermal conductivity is pressure- and crystallinity-dependent. Specific heat depends on temperature and crystallinity. Latent heat is released according to the crystallization progress and temperature. Deviatoric elastic stress is evolved in the solidified material. The prediction of the cavity pressure evolution is used for the assessment of the solution quality because it is experimentally available and governs the residual stress development. Insight into the thermomechanical conditions is provided with through-thickness plots of pressure, temperature and cooling rate at different levels of crystallinity. The code and simulation setup are made openly available to further the research on the topic.


Introduction
Injection molding simulation research requires a considerably advanced computer code, which is particularly true when including crystallization modeling. Many publications made use of in-house codes that were not made public, while the industrial injection molding solvers specialize for the industrial needs and prove to be restrictive for research of the thermomechanical phenomena in injection molding. The commercial general purpose computational fluid dynamics (CFD) codes may allow customization for conducting an injection molding simulation, but the undertaking is challenging and likewise results an in-house code.
Injection molding simulation has been the subject of ongoing research for decades. Kennedy and Zheng [1] have thoroughly reviewed the history of injection molding simulation publications. Based on the aim of predicting part geometry, the research was focused on the residual stress prediction. An important publication was contributed by Baaijens [2] who developed a model capable of describing the filling and packing stages as well as the final residual stresses. They highlighted the effect of mold compliance on the cavity pressure evolution. Later, the constitutive modeling was advanced with an advanced viscoelastic material model by Chang and Chiou [3], a three-dimensional finite volume method was applied [4] and residual stresses were analyzed both experimentally and numerically for amorphous and crystalline polymers [5]-Guevara-Morales and Figueroa-López [6] published a thorough review of the research on the residual stresses.
An important contribution to injection molding research was made by the development of the UNISA code, described by Pantani et al. [7], which scientifically tackled the crystallization specifics in addition to the general problem of injection molding. The code is advanced in terms of physics, but does not appear to be publicly accessible. It was used by De Santis et al. [8] for predicting the shrinkage of injection molded isotactic polypropylene (iPP) and inspecting the role of holding time and pressure. Pantani et al. [9,10] also investigated crystallization modeling in injection molding simulation, compiling a comprehensive review of the subject. The polymer of interest was iPP. Zheng et al. [11] modeled flow induced crystallization while also accounting for the colorant content, likewise investigating iPP. Zheng et al. [12] reviewed the modeling aspect in a dedicated chapter while Janeschitz-Kriegl [13] compiled a thorough introduction to the physics of polymer crystallization and provided material data for different polymers, demonstrating the complexity of the phenomenon form the experimental and numerical perspective.
An application-oriented study of high-density polyethylene (HDPE) injection molding was performed by Kabanemi et al. [14]. The study focused on the part geometry prediction in the scope of solid mechanics, while the fluid thermo-mechanics was of secondary importance. Kamal et al. [5] simulated injection molding of HDPE using a finite volume code McKam4. The software was capable of describing crystallization, but has not been referenced in later studies and does not appear to be publicly available. Ilinca and Hetú [15] simulated gas-assisted injection molding of HDPE where they mostly elaborated on the numerical procedure.
While the studies concerned with crystallization in injection molding simulation focused on iPP, HDPE has received far less attention. The published work tended to employ commercial codes or inhouse codes. On this basis, we formulate two aims of this study: To overcome the hindrance of commercial code lock-in and aid the research community cooperation, we have customized the open-source general purpose CFD code OpenFOAM [16] and made it publicly available as openInjMoldSim [17]. The solver is a modification of compressibleInterFoam of OpenFOAM v3.0.1 [16]. Within the scope of our previous work [18], we demonstrated this approach on modeling volumetric relaxation of amorphous polystyrene during injection molding using a modified openInjMoldSim code, named openInjMoldDyMSimAmClr [19]. The present publication demonstrates a code version for simulating injection molding of HDPE with crystallization, denoted as open-InjMoldDyMSimCr [20]. This highlights the research potential of an open-source injection molding code that can be tailored to the specific research needs.
The numerical solver is suitable for non-simplified three-dimensional geometry and incorporates three influences on the cavity pressure evolution: mold compliance, pressuredependent thermal contact with the mold and the crystallization dependent specific volume. This allows comparing the prediction to the experimental pressure evolution.
The course of this work (section: Materials and methods) consists of three major parts. The experimental results are outlined and the reader is referred to the publication with further details. The crystallization model introduces the theoretical background of the HDPE crystallization model. Subsequently, we describe the formulation of a numerical model of injection molding to assess the quality of the prediction of the pressure evolution by comparing it to the experimental data. The results are reviewed and discussed, offering quantitative insight into the thermo-mechanics of an injection molded HDPE product. The paper is concluded with the implications for the research field and directions for future work.

Experimental Investigation
The experimental counterpart of the examined case consists of pressure evolutions that were measured and published before [21]. The temperature of the mold cooling water was set to 50 • C, leading to about 52 • C at the surface of the cavity while the melt temperature was set to 220 • C. The crystalline plastic HDPE Sabic M80064 was used. The experimental mold cavity is shown in Figure 1 with the pressure measurement positions denoted by P0, P1, P3 and P4. The pressure measuring method was indirect, i.e., using rods installed in the mold that transmitted the force from the cavity to the sensor. The relevant pressure evolutions are reported together with the numerical results later on.

HDPE Thermal and Volumetric Behavior
Crystallization influences the thermal and mechanical conditions. The transition from melt to solid implies a change in the specific heat, a release of latent heat [22] and a decrease in the specific volume. An increase in thermal conductivity due to solidification has also been measured [23], and is taken into account here.

Crystallization
The crystallinity degree was evolved using the Kolmogofoff-Avrami-Evans model [1] expressed in the form of the Scheinder's differential equation system [24]. The fictive crystalline volume field ϕ is integrated via intermediate variables ϕ 1 and ϕ 2 according to, where N k is the nuclei density and G k the spherulite lineal growth rate. The volume fraction of the solidified material is calculated as: The mass fraction of the solidified material is taken as the relative crystallinity, where ρ (m) and ρ (s) are the melt and solid density, respectively. The pressure influence modeling was adopted from Zuidema et al. [25], with a decrease in the effective temperature as, where b 6 is adopted from the Tait equation while T c is the temperature input to the crystallization model. The nucleation density (nuclei number per volume unit) description was adopted from Koscher and Fulchiron [26], where T 0 m is the HDPE melting temperature at 419 K [27]. The free parameters N 0 and ∆T 0 were identified in the scope of this work by matching the model crystallization rate with the HDPE calorimetry measurements [28,29] (Table 1).
The spherulite lineal growth rate was described according to Mandelkern et al. [30] as reported by Van Krevelen and Nijenhuis [27] to which an additional term was added as suggested by Mandelkern [31], with R p as the universal gas constant. The parameters G k,0 , C 3 and E D , were provided for polyethylene by Van Krevelen and Nijenhuis [27] ( Table 1). The parameter r c phenomenologically accounts for the fractions of untransformable material in the kinetics [31] and was identified in this work for HDPE by matching the calorimetry data to the crystallization model output.
The temperature at half-crystallization T cr = T(ξ = 0.5)0 was plotted against available experimental data for a broader range of cooling rates ( Figure 2) and compared to available fast cooling rate measurements on polyethylene [28,29,32,33]. Zhuravlev et al. [29] tested an HDPE with a molecular mass of 500 kDa, while Androsch et al. [33] tested ultra-high molecular weight polyethylene (UHMWPE) and obtained similar half crystallization times. Toda et al. [28] examined a polyethylene of 50 kDa.

Specific Volume
The specific volume is expressed as, where v (s) and v (m) are the solid and melt specific volume, respectively [12,34]. The v (s) and v (m) are identifiable from the two-domain Tait equation (see Zheng et al. [12]) which is also derived based on the rule of mixture, and The subscripts "m" and "s" denote the b parameters' values for the melt and solid domains, respectively, with the universal constant C = 0.0894. The parameters are listed in Table 2 for HDPE Sabic M80064 [35].
Integrating the crystallization degree at constant cooling (Equations (1) and (3)), the specific volume of HDPE (Equation (7)) is plotted against temperature for different cooling rates in Figure 3, revealing the crystallization effect. As the temperature is considered independent, the thermal properties do not affect this result.
The subscripts "m" and "s" denote the parameters' values for the melt and solid domains, respectively, with the universal constant = 0.0894. The parameters are listed in Table 2 for HDPE Sabic M80064 [35]. Integrating the crystallization degree at constant cooling (Equations (1) and (3)), the specific volume of HDPE (Equation (7)) is plotted against temperature for different cooling rates in Figure 3, revealing the crystallization effect. As the temperature is considered independent, the thermal properties do not affect this result.

Specific Heat
Latent heat is released with crystallization and has a strong influence on the thermal conditions. Gaur and Wunderlich [22], described the apparent specific heat by lumping the latent into the to the mixture equation,

Specific Heat
Latent heat is released with crystallization and has a strong influence on the thermal conditions. Gaur and Wunderlich [22], described the apparent specific heat by lumping the latent into the to the mixture equation, where c p is the specific heat of the mixture at constant pressure, c c p and c a p are the specific heat at constant pressure for the crystalline, and the amorphous phase, respectively, w c is the mass fraction of the crystalline phase, and ∆H f is the heat of fusion. In this work, the crystalline mass fraction is calculated according to, by assuming the semi-crystalline phase to contain a fraction of w c ∞ = 0.66 of the crystallized material. In the numerical implementation of Equation (10), the first two terms are used to describe the specific heat with the last term transferred to a heat source adopting the values recommended by Gaur and Wunderlich [22]. The resulting apparent specific heat is shown in Figure 4 for different cooling rates.

Injection Molding Model Formulation
The development of crystallization was investigated in a simulation of the injection molding filling and packing stages. This is a problem of a laminar, non-isothermal, compressible flow of two immiscible fluids according to the volume of fluid method. The fluid mechanics problem is based on the conservation equations complemented by the constitutive equations. A system of partial differential equations is obtained and discussed from the geometrical aspect. The code and the numerical setup are publicly available under the name openInjMoldDyMSimCr [20] as a further modified version of openInjMoldSim [17].

Conservation Equations
The code employs the finite volume method to solve the differential conservation equations for the mass, momentum and energy. The conservation of mass relates the density ρ to the velocity field u k , ∂ρ ∂t where Einstein's notation is used with k = 1, 2, 3 for the Cartesian coordinates x k . The conservation of momentum relates the velocity field to the stresses σ ij , where the notation for the material derivative is used: The stresses are composed of the hydrostatic pressure p, the elastic contribution τ e ij in the solidified material and the viscous contribution τ v ij as, with the Kronecker delta tensor δ ij . Finally, the energy conservation is imposed with the equation, with the specific heat c p , thermal conductivity k and the rate of strain tensor: The thermal conductivity was adopted from Dawson et al. [23] as shown in Figure 5. The depicted thermal dependence only approximately illustrates the effect of crystallization. Note that Dawson et al. [23] performed the measurements during heating of HDPE which lead to a phase change at higher temperatures than during cooling. Figure 5. The thermal conductivity as measured [23] and modeled (thick horizontal lines).

Viscosity
The frequently used Cross-WLF melt viscosity model was employed, . γ τ * n−1 (19) with the model parameters listed in Table 3 as obtained from [35] and corrected using rotational rheometry.

Constitutive Modeling
The viscosity and specific volume were modeled as described in the previous section with the viscous deviatoric stress was calculated as, (20) where D d ij is the deviatoric component of the rate of deformation tensor. The melt was assumed to solidify when reaching the maximum viscosity of 0.5 MPas, at which point the deviatoric elasticity was onset according to the Upper Convected Maxwell model [36] with an infinite relaxation time, where τ e ij is the elastic deviatoric stress in Einstein's notation, u i is the velocity vector, x i are the Cartesian coordinates, G = 200 MPa is the shear modulus [37] and D d ij is the deviatoric rate of deformation tensor.
The growing fraction of the crystalline content in the melt gradually leads to flow cessation. The increase in viscosity is a challenging phenomenon to investigate experimentally. The viscosity was scaled due to crystallinity using the empirical factor, as suggested by Titomanlio et al. [38] (Table 4).

Geometry and Boundary Conditions
The cavity was described by the longitudinal cross-section displayed in Figure 6. The measured P1 pressure was used as the boundary condition. This eliminated the gate from the computational domain where the flow is not planar. The air was released through the outlet during the filling stage which was closed during the packing stage. The melt filling temperature was 220 • C and the mold temperature 50 • C. The heat transfer coefficient h on the mold walls was modeled as pressure dependent with a linear interpolation between the points (p, h) = 0 MPa, 1020 W/ m 2 K and (100 MPa, 5384 W/(m 2 K)) as identified from the work of Delaunay et al. [39]. The compliance of the cavity was set to c = 0.75 µm/MPa, similarly as in [40]. The finite volume mesh consisted of 64 cells in the thickness direction and 2400 cells in the length direction.

Results and Discussion
The calculated pressures at P2 and P3 positions are compared to the measured values ( Figure 7). The P1 pressure was imposed at the inlet according to the measurements. This result offers insight into the quality of the pressure gradient and filling rate prediction. The evolution of the packing pressure at these positions is depicted in Figure 8. The solidification time appears to have been well matched with the experimental values, while pressure gradient during post-filling was found to be slightly over-predicted. Further investigation of the post-filling phase pressure gradient would involve an experimental and numerical investigation of the D 3 parameter in the Cross-WLF Equation (18) as it increases the high-pressure viscosity. At the start of the cooling a temperature increase above the initial melt temperature is predicted due to adiabatic compression. The latent heat release was most evident at 120 • C in the mid-thickness as shown in Figure 9 for the experimental pressure positions. A plateau developed as typically in quenching of crystalline materials. A slight temperature increase was predicted at the start of the plateau which could be realistic [41]. A far faster temperature drop was found at the mold contact ( Figure 10) where a temperature rise during packing developed as a consequence of the decreasing heat transfer coefficient at dropping pressure (see [39,42]). A far faster temperature drop was found at the mold contact ( Figure 10) where a temperature rise during packing developed as a consequence of the decreasing heat transfer coefficient at dropping pressure (see [39,42]). At 4 s the crystallization on the P2 position was almost complete (Figure 11) and the graphs indicate that the crystallization near the mold surfaces was far more rapid.  At 4 s the crystallization on the P2 position was almost complete ( Figure 11) and the graphs indicate that the crystallization near the mold surfaces was far more rapid. The times of reaching different relative crystallinity degrees are depicted in Figure 12, indicating that the center of the cavity took more than a second to advance the relative crystallinity from 1% to 50%, unlike the surface layers with rapid crystallization. Temperature at solidification is an important piece of information when modeling is simplified by assuming the so called "no flow temperature" [43]. Figure 13 depicts the temperature through the thickness when the 50% relative crystallinity was reached. It was found to be confined between 110 • C and 120 • C.
The time derivative of temperature at different values of relative crystallinity ( Figure 14) reveals that the near surface layers were cooling at rates up to −130 K/s, while the latent heat lead to a brief and rapid increase in temperature. Temperature at solidification is an important piece of information when modeling is simplified by assuming the so called "no flow temperature" [43]. Figure 13 depicts the temperature through the thickness when the 50% relative crystallinity was reached. It was found to be confined between 110 °C and 120 °C. The time derivative of temperature at different values of relative crystallinity ( Figure 14) reveals that the near surface layers were cooling at rates up to −130 K/s, while the latent heat lead to a brief and rapid increase in temperature.  The pressure change during crystallization was smaller in the near-mold material ( Figure 15). This result provides insight into the impact of assuming instant solidification in calculating residual stresses.

Conclusions
The use of an open-source code for injection molding simulation with crystallization modeling was demonstrated. A validation was provided with the comparison of the predicted and the measured cavity pressure evolutions. While other research focused on iPP, the HDPE crystallization was inspected. The Kolmogorof-Avrami-Evans crystallization model was calibrated according to fast scanning calorimetry studies and insight into relevant cooling rates and solidification temperatures was gained. This can, for instance, serve as guidance in selecting the no-flow temperature or adjusting the two-domain Tait equation in industrial simulation.
The injection molding simulation research can be conducted using OpenFOAM, an open source CFD package, allowing the research to focus to the physics by utilizing the established CFD workflows. The amount of the required coding to develop similar solutions was thus greatly reduced.
The use of OpenFOAM has important implications for further research. Being a state-of-the-art general purpose CFD code while also open source it allows for advanced numerical modeling, with options for tailoring the numerical solution by altering the linear system of equations' solvers and appropriate differencing schemes. Advanced modeling approaches of dynamic meshing allowed modeling a deformable cavity but could be used to model compression molding etc.
The utilized code and simulation setup are publicly available including the code for generating the graphical results. This allows the researchers to probe the influence of any of the listed parameters, including cavity thickness or flow-path length. The modeling can be further developed. The crystallization model can be updated to include flow induced crystallization and to predict non-spherulitic morphology. As for the general features, modeling of gates, venting and in-mold shrinkage would contribute a great deal to its industrial applicability. Ultimately, it could be computationally optimized to also serve the industry where accuracy is often sacrificed to reduce computational times.