Multiscale Numerical Study of 3D Polymer Crystallization during Cooling Stage

Abstract

We aim to study the behavior of polymer crystallization during cooling stage in injection molding more accurately, the multiscale model and multiscale algorithm proposed in our previous work (Ruan et al., 2012) have been extended to the 3D polymer crystallization case. Our multiscale model incorporates two distinct length scales: a coarse grid for the heat diffusion and a fine grid for the crystal morphology evolution (nucleation, growth, and impingement). Our multiscale algorithm couples the different methods on different length scales, namely, the finite volume method (FVM) on the coarse grid and the pixel coloring method on the fine grid. By using these multiscale model and multiscale algorithm, simulations for 3D polymer crystallization are carried out. Macroscopic variables, for example, temperature, relative crystallinity, as well as the microscopic structural characters, for example, crystal morphology development, and mean size of spherulites, are investigated at various cooling conditions. We also show the importance of coupling heat transfer with crystallization as well as 3D numerical studies.

1. Introduction

Nowadays, semicrystalline polymers play an important role in industry due to their advantages of enhanced mechanical properties, ease of manufacturing, and so forth [1]. These polymers are typically processed by injection molding processing techniques to form the products. As one of the most widely used polymer processing techniques, injection molding mainly consists of filling, packing, and cooling [2]. Among these, cooling stage is the most significant part not only because it accounts for the largest part of the total injection molding cycle, but also for the crystal morphology formed during crystallization which is known as a crucial factor to the physical and mechanical properties of products. The crystal morphology also known as the microstructure is in turn determined by a number of processing conditions during the cooling stage. Thus, it is of great importance to accurately model the crystallization during cooling stage and predict the final crystal morphology formed under different processing conditions [1].

Polymer crystallization is a typical multiscale processing: the molecular chains fold together and form ordered regions called lamellae, which compose larger spheroidal structures named spherulites [1]. Since the polymer crystallization has a span of about 10¹⁰ (from the molecular chain length 10⁻⁸∼10⁻¹⁰ m to the product length 10⁻²∼10⁰ m [1]) in length scale, it is hard to carry out a full-scale simulation to bridge single molecular chain folding to the final product. Therefore, to model the polymer crystallization, one should first choose a suitable length scale. Spherulite level (10⁻³∼10⁻⁵ m [1]) is proper because it is the single largest feature of polymer microstructures, in addition, it is relatively easier to link the heat transfer in the macroscopic level (10⁻²∼10⁰ m [1]).

To date, a number of numerical investigations have appeared on spherulitic crystallization in polymer melts during cooling. Charbon and Swaminarayan [3] proposed a multiscale algorithm for the simulation of polymer crystallization. They incorporated front-tracking techniques with nucleation, spherulite impingement, and latent heat release as well as heat diffusion to predict the evolution of microstructures and relative crystallinity during cooling stage. Huang and Kamal [4] presented a physical model for polymer crystallization. They coupled a nonconserved envelop vector field of the local macroscopic growth domains with a vector field of the local mesoscopic lamellae directors. Prabhu et al. [5] proposed a coupling finite-element method and the cell model to the numerical study of polymer crystallization. The importance of micro-macro-coupling was analyzed and explained. Recently, Ruan et al. [6] proposed a coupling FVM and the pixel coloring method to deal with crystallization in short-fiber-reinforced composites. Crystal morphology evolution was explicitly shown by the pixel coloring method. However, it should be pointed out that the aforementioned papers only deal with the 2D case which is a qualitative simplification for the actual problem. Moreover, works of Charbon and Swaminarayan [3] as well as Huang and Kamal [4] only deal with the imposed temperature filed, the effect of processing conditions has not been investigated.

In the numerical study of 3D polymer crystallization, Raabe [7], Raabe and Godara [8], and Lin et al. [9] were the pioneers. The cellular automaton method was constructed to deal with the detailed development of crystal morphology. However, their works are only concentrated on the simple isothermal case. So far, no literature has been found to deal with the crystal morphology development with the heat transfer in 3D numerical study.

The objective of this article is to extend the multiscale model and the multiscale computational method proposed in our pervious work [6] to the 3D numerical study of polymer crystallization during cooling. Since the injection molding is 3D in realistic situation, 3D numerical study can predict the temperature distribution and the crystallization behavior more accurately. Our multiscale model incorporates two distinct length scales to simulate the crystallization: a coarse grid for the heat diffusion and a fine grid for capturing the crystal morphology formation. FVM is employed on the coarse grid to solve the energy equation while the pixel coloring method is applied on the fine grid to track the crystal morphology development. The present paper is organized as follows: in Section 2, the multiscale model and the coupled computational method are introduced; in Section 3, results and discussion are given; finally the conclusions are drawn in Section 4.

2. Multiscale Model and Multiscale Algorithm

During the cooling stage, the process of crystallization is coupled with the heat transfer which makes the modeling more difficult and complex. On the one hand, it is well known that the kinetic parameters of nucleation and growth in microscopic scale are strongly dependent on the temperature or processing conditions. On the other hand, the crystallization is an exothermic process which releases the heat and affects the thermal field in macroscopic level.

2.1. Thermal Field in the Macroscopic Level

Thermomechanical histories are of great importance in the estimate of the final product properties. The corresponding simplified energy equation is [6] 𝜌𝑐_𝑝𝜕𝑇𝜕𝑡=∇⋅𝑘^𝑝∇𝑇+𝜌Δ𝐻𝜕𝛼𝜕𝑡,(2.1) where 𝜌 is the material density, 𝑐_𝑝 is the thermal capacity, 𝑇 is the temperature field, 𝑡 is the time, 𝑘_𝑝 is the thermal conductivity, Δ𝐻 is the total enthalpy released during crystallization, and 𝛼 is the relative crystallinity. The last source term is the contribution of the latent heat released by the crystallization, which plays the role in macro-micro-coupling.

In the accurate modeling, it is important to consider the thermal properties as a function of temperature and relative crystallinity. The thermal capacity, thermal conductivity can be described by the “mixing rule” of the solid state and the liquid state values to get [10] 𝑘𝑐_𝑝𝑝=𝛼𝑐_psc+(1−𝛼)𝑐_pa,=𝛼𝑘_psc+(1−𝛼)𝑘_pa,(2.2) where 𝑐_psc and 𝑐_pa are the thermal capacity of semicrystalline phase and amorphous phase, respectively, 𝑘_psc and 𝑘_pa are the thermal conductivity of semicrystalline phase, and amorphous phase, respectively. Also 𝑐_psc, 𝑐_pa, 𝑘_psc, and 𝑘_pa are temperature dependent variables, which can be modeled as a linear temperature dependency functions [10].

2.2. Crystal Morphology Evolution in the Microscopic Level [6]

Crystallization is a mechanism of phase change in semicrystalline polymeric materials [11]. It represents a mix between nucleation, growth, and impingement. When the temperature of semicrystalline polymeric melt is decreased below its crystallization temperature, “nuclei” appears randomly in space and time; then the appeared nuclei start to grow to form what is referred to as “spherulites.” As the time progresses, the spherulites grow until they hit another crystal and stop growing which is called “impingement.” If all the spherulites are impinged with other spherulites, the crystallization process is finished.

2.2.1. Nucleation

Polymer nucleation is an important factor which affects the final morphology. Usually, the nucleation may vary material from material. In the numerical study, one may adopt an empirical nucleation relation which is derived from fitting the experiment data. Here, we use the following relation of nucleation density [12, 13]: 𝑁(𝑇)=𝑁₀exp[𝜑Δ𝑇],(2.3) where Δ𝑇=𝑇^0𝑚−𝑇 is the supercooling temperature with 𝑇^0𝑚 being the equilibrium melting temperature, 𝑁₀ and 𝜑 are the empirical parameters. This nucleation density relation was also used in our pervious work [6] for the description of nucleation density in polymer bulk.

2.2.2. Growth and Impingement

Growth rate is another important factor which affects the development of morphology. Here, we adopt the Hoffman-Lauritzen theory [14] to describe the spherulite growth rate: 𝐺(𝑇)=𝐺₀exp−𝑈^∗𝑅_𝑔𝑇−𝑇_∞exp−𝑇Δ𝐾𝑇^𝑔𝑓,(2.4) where 𝐺₀ and 𝐾_𝑔 are constants, 𝑈^∗ is the activation energy of motion, 𝑅_𝑔 is the gas constant, 𝑇_∞=𝑇_𝑔−30 with 𝑇_𝑔 is the glass transition temperature, and 𝑓=2𝑇/(𝑇^0𝑚+𝑇).

With the growth of spherulites, it is inevitable to meet with the “impingement.” Impingement is happed in spherulites which contact with their neighbors or the walls. Different spherulites impinge to form grain boundaries. With the help of grain boundaries, it is possible to calculate the mean size of spherulites.

2.3. Macro-Micro-Algorithm

In our paper, we assume that the temperature field is in the macroscopic level while the crystal morphology evolution is in the microscopic level (see Figure 1). Therefore, it is important to build up a macro-micro-algorithm. Here we extend the algorithm of coupling of FVM [15] and the pixel coloring method [6, 11, 16] proposed in our pervious work [6] to the 3D numerical study.

Figure 1 

Schematic representation of different length scale for computation.

In our algorithm, the domain is first divided by a coarse grid. FVM with cell vertex scheme is used on this coarse grid to calculate the temperature. Then, each coarse grid is subdivided into a number of cubes to form the fine grid. The pixel coloring method is employed on this fine grid to track the development of crystal morphology. We assume the temperature on each coarse grid is the same which determines the nucleation and growth rate of spherulites through (2.3) and (2.4). On the other side, with the evolution of crystallization (nucleation, growth, and impingement) on the fine grid, the relative crystallinity changes which affects the temperature through (2.1). This realizes the macro-micro-coupling.

On the coarse grid, FVM is used to solve the energy equation. The reason why we use FVM is because this method uses the control volume which is more like the cell in the statistics for the microscopic information. Figure 2(a) shows a 3D control volume for FVM. Here, 𝑃 is the computational point with the neighbor points 𝐸, 𝑊, 𝑁, 𝑆, 𝑇, 𝐵 and the interface points 𝑒, 𝑤, 𝑛, 𝑠, 𝑡, 𝑏. See Figure 2(b) which gives a 2D gird. We use the first-order forward-time and second-order central-space scheme to discrete the energy equation (2.1) to get 𝜌𝑐_𝑝𝑇^𝑃𝑛+1−𝑇^𝑛𝑃Δ𝑡=𝑘^𝑝𝑇^𝑛𝐸−2𝑇^𝑛𝑃+𝑇^𝑛𝑊Δ𝑥²+𝑇^𝑛𝑁−2𝑇^𝑛𝑃+𝑇^𝑛𝑆Δ𝑦²+𝑇^𝑛𝑇−2𝑇^𝑛𝑃+𝑇^𝑛𝐵Δ𝑧²+𝜌Δ𝐻𝛼^𝑛𝑃−𝛼^𝑃𝑛−1Δ𝑡.(2.5) The choice for such a scheme is because of numerical simplicity. In our algorithm, we first calculate the temperature then calculate the relative crystallinity, therefore, the relative crystallinity is obtained with a delay that is why we use 𝜌Δ𝐻(𝛼^𝑛𝑃−𝛼^𝑃𝑛−1)/Δ𝑡 instead of 𝜌Δ𝐻(𝛼^𝑃𝑛+1−𝛼^𝑛𝑃)/Δ𝑡 for the last term in (2.5). It should be noted that this energy equation calculation is to obtain temperature which is the input of nucleation and growth of spherulites on fine grid.

(a)

(b)

Figure 2 

Schematic of a control volume (a) 3D grid (b) 2D grid (profile of 3D grid).

On the fine grid, the pixel-coloring method [6, 11, 16] is implemented to capture the growth front of a population of spherulites. Here, we will not present the detailed implementation of the pixel-coloring method and refer our previous work [6] for more details. The statistical character of the evolution of crystal morphology on fine grid is the relative crystallinity which is the input of temperature calculation on coarse gird. In our study, the relative crystallinity can be calculated by 𝛼=numberofcellsthathavebeentransformedtotalnumberofcells.(2.6) Since the algorithm on fine grid can explicitly show the details of spherulites, here, we define another important factor, the mean radius of each spherulite, as 𝑅=³3𝑉(4𝜋)(2.7) with 𝑉 the volume of the spherulite which can be calculated with the number of cells occupied by the spherulite and the cell size. It should be mentioned that the mean radius of spherulite directly affects the stress of the final product through the Hall-Petch relation [11].

Figure 3 shows the macro-micro-algorithm used for our simulation.

Figure 3 

Multiscale algorithm.

3. Results and Discussion

3.1. Problem Formulation

3D polymer crystallization during cooling stage is studied here. Figure 4 describes the computational geometry. The length of the mold cavity is 8 mm while its width and height are 4 mm respectively. We assume the walls 𝑦=0 mm, and 𝑧=0 mm experiencing the constant cooling rate operation and the temperature boundary conditions are set to be 𝑇=𝑇₀−𝑐×𝑡, with 𝑇₀ the initial temperature and 𝑐 the cooling rate. The other boundaries are assumed adiabatic and the temperature boundary conditions are set to be 𝜕𝑇/𝜕𝐧=0 with 𝐧 the outward unit normal vector. It should be mentioned that this geometry it related to the think-wall parts in injection molding which is because the ratio of width/height of the cross-section is less than 10 [2].

Figure 4 

Computational geometry.

The material we considered here is iPP and the parameters used in the simulation are [6, 12, 13]: 𝑁₀=17.4×10⁶/m³, 𝜑=0.155, 𝐺₀=2.1×10¹⁰𝜇m/s, 𝑈^∗/𝑅_𝑔=755K, 𝐾_𝑔=534858K², 𝑇^0𝑚=467K, 𝑇_𝑔=266K, 𝜌=900kg/m³, 𝑐_𝑝=2.14×10³J/kg/K, 𝑘_𝑝=0.193W/m/K, Δ𝐻=107×10³J/kg. Unless otherwise stated, the other parameters are set to be 𝑇₀=470K, 𝑐=2K/min. Due to the lack of material parameters, here we neglect the thermal difference between the solid state and the liquid state.

We test three meshes in this problem computation, namely, Mesh1: 8×4×4 for coarse grid and 100×100×100 for fine grid; Mesh2: 12×6×6 for coarse grid and 100×100×100 for fine grid; Mesh3: 8×4×4 for coarse grid and 150×150×150 for fine grid. Results of temperature and relative crystallinity at the intersection line of 𝑥=4 mm plane and 𝑧=4 mm plane (line: 𝐴𝐵, see Figure 4) are compared in Figure 5. As we can see that the result obtained on Mesh1 agrees with the results obtained on Mesh2 and Mesh3 very well. As we know, the refinement of either coarse grid (Mesh2) or fine grid (Mesh3) can lead to more accurate of the final results, it brings in a large of storage which affects the efficiency and challenges the computer. When we balance the accuracy with the efficiency, we find that Mesh1 is proper. Therefore, in our following study, Mesh1 is used and its coarse grid size is 1mm×1mm×1mm and the fine grid size is 10𝜇m×10𝜇m×10𝜇m. Here, we will not explain the phenomenon of temperature and relative crystallinity obtained in Figure 5 and describe it in the next section.

(a)

(b)

Figure 5 

Evolution of temperature and relative crystallinity at the intersection line of 𝑥=4 mm plane and 𝑧=4 mm plane (line: 𝐴𝐵).

3.2. Importance of Macro-Micro-Coupling

3.2.1. Temperature and Relative Crystallinity Evolution in Macroscopic Level

Evolution of temperature and relative crystallinity during cooling at the plane of 𝑥=4 mm are shown in Figure 6(a). As we can see in the figure of temperature evolution, a quasi-isothermal plateau is formed at the positions near the adiabatic boundary (𝑦=4 mm, 𝑧=4 mm). This plateau is caused by the latent heat released by crystallization. We also reported this plateau in our pervious work [6] which deals with the 2D short-fiber-reinforced system. The figure of relative crystallinity evolution tells us that the polymer in the skin layer crystallizes much earlier and faster than that in the core layer. This can be explained that the supercooling temperature in the skin layer is much higher than that in the core layer due to the temperature difference which is the fundamental reason for the crystallization.

(a)

(b)

Figure 6 

Evolution of temperature and relative crystallinity (a) considering the latent heat (b) without considering the latent heat.

Latent heat released by the crystallization is the bridge in macro-micro-coupling. To determine the importance of the contribution of the latent heat, the results of our simulation are compared with the results of model which ignores the latent heat. Figure 6(b) shows the results of model without consider the latent heat. We can see that there is a large difference between these two models especially at the positions near the adiabatic boundaries. Moreover, the results of model without considering latent heat do not develop the quasi-isothermal plateau at the positions near the adiabatic boundary. This is the error caused by ignoring the latent heat. Since the nucleation and growth rate of the spherulites is completely determined by the temperature, therefore, to predict crystal morphology more accurately, the model of temperature is necessary to include the influences of the latent heat.

3.2.2. Crystal Morphology in Microscopic Level

Evolution of crystal morphology at the control volume of 𝑥=4 mm is shown in Figure 7. We use the fronts of spherulitics in surface (or slices) to show the crystal morphology evolution. Here, the white region is the polymer melt and the colored region is the spherulitics. Different spherulitics are distinguished by different colors. Figure 7 clearly shows that crystals first appear in the skin layer, with the time evolutes, they gradually appear to the core layer. This is consistent with the change tendency of relative crystallinity (see Figure 6(a)) which is a macrostatistical character of crystal morphology evolution in microscopic level. The results presented here show qualitatively agreement with our pervious work [6] for the 2D short fiber reinforced polymer solidification.

(a)

(b)

(c)

Figure 7 

Evolution of crystal morphology at the control volume of 𝑥=4 mm (a) 𝑡=2360 s, (b) 𝑡=2400 s, (c) 𝑡=2600 s.

Final morphology in the skin layer and in the core layer is compared in Figure 8. We choose skin layer as a representative control volume of point with the coordinate (4 mm, 1 mm, 1 mm) while core layer as a representative control volume of point with the coordinate (4 mm, 3 mm, 3 mm) which is illustrated in Figure 4. It is not so obvious for us to distinguish the skin morphology from the core morphology. However, the size of spherulitics in the core layer is somewhat bigger than that in the skin layer.

(a)

(b)

Figure 8 

Crystal morphology at different positions (a) skin (b) core.

Mean radius of each spherulite and the distribution of spherulite size at different positions are shown in Figure 9. We should mention that we do not consider the “phantom nuclei” in the comparison. Figure 9 shows that spherulite size is relative smaller in the position which is close to skin layer. Since the spherulites size directly affects the mechanical properties of the products according to Hall-Petch relation [11], it can be concluded that mechanical properties vary from skin to core.

(a)

(b)

Figure 9 

Spherulite size and distribution of spherulite size at different positions.

3.3. Importance of 3D Simulation

In order to highlight the importance of 3D simulation, we also give the comparison of the results obtained in 2D case with the 3D case. Since the problem we studied is a quasi-three dimensional problem, it can be also reduced as a two-dimensional one. Here, we consider the profile of 𝑥=4 mm (see Figure 4) as our 2D studied cavity with the boundaries 𝑦=0 mm and 𝑧=0 mm experiencing the cooling rate operation and the other boundaries been assumed to be adiabatic. We will mention that the only difference between 3D case and 2D case is that the heat transfer in 𝑥 direction is considered in 3D case while its effect is ignored in 2D case. In our study, the nucleation density in 2D and 3D is related with the following statistical formulation [3]: 𝑁_2D=1.458×(𝑁_3D)^2/3 and the growth rate of spherulites in 2D and 3D are assumed to be the same.

Figure 10 shows the evolution of temperature and relative crystallinity which are obtained in 2D simulation. Through comparing with Figure 6(a), it is observed that the temperature away from the cooling boundaries is higher in 2D case than 3D case at the same time; particularly, the maximum temperature divergence is about 1 K in the core layer. Moreover, the internal relative crystallinity is also higher in 2D case than 3D case; in particular, the maximum relative crystallinity divergence is about 0.15 in the core layer.

(a)

(b)

Figure 10 

Evolution of temperature and relative crystallinity (2D results).

Evolution of crystal morphology at the whole computational geometry in 2D case is shown in Figure 11. Compared with Figure 7, we can find that the tendency of microstructure formation in 2D case is the same as 3D case. To emphasize the difference between 2D case and 3D case, we also investigate the spherulite size and its distribution at different positions of the cavity. For the sake of simplicity, we also choose the “skin” and “core” control volume as our comparison cells as illustrated in Figure 4. Figure 12 shows the spherulite size and its distribution as obtained in 2D simulation. By the comparison with Figure 9, we find that the trends of spherulite size distribution in skin and core are identical in 2D case and 3D case, however, in 2D case, the distribution of spherulite size is more concentrated.

(a)

(b)

(c)

Figure 11 

Evolution of crystal morphology in 2D case (a) 𝑡=2360 s (b) 𝑡=2400 s (c) 𝑡=2600 s.

(a)

(b)

Figure 12 

Spherulite size and distribution of spherulite size at different positions (2D results).

The comparison in this section tells us that the macroscopic and microscopic values obtained in 2D simulation are quantitatively different from that in 3D simulation. Therefore, if the models, algorithms, and computational conditions are allowed, we should consider the 3D simulation in order to obtain the more precise results.

3.4. Effects of Thermal Condition on the Crystallization

Effects of cooling rate and initial temperature are also investigated in this 3D polymer crystallization case. Without loss of generality, we here choose the “core” control volume as our cell for showing the results.

3.4.1. Effects of Cooling Rate

We change the boundary conditions by varying the cooling rate as 1 K/min, 2 K/min, 5 K/min in order to study the effects of cooling rate.

Figure 13 shows the relative crystallinity evolution and the distribution of spherulite size in the core layer for different cooling rate. According to Figure 13, the higher cooling rate causes an acceleration of the crystallization and a reduction of spherulite size. We will mention that the effects of cooling rate are similar as our pervious work [6] where the 2D short-fiber-reinforced system are studied. In addition, this tendency is also reported in the experimental work by Zheng et al. [17].

(a)

(b)

Figure 13 

Evolution of relative crystallinity and the distribution of spherulite size in core layer for different cooling rates.

Thus, we will conclude that if the designer wants to obtain the smaller size of the spherulites, he shall impose a relatively higher cooling rate boundary condition.

3.4.2. Effects of Initial Temperature

We varying the initial temperature as 470 K, 480 K, 490 K in order to study the effects of initial temperature.

Figure 14 shows the relative crystallinity evolution and the distribution of spherulite size in the core layer with different initial temperature. It is observed that the higher initial temperature leads to a deceleration of crystallization and minor effect on the crystal morphology. This is also in agreement with our pervious work [6] in 2D case.

(a)

(b)

Figure 14 

Evolution of relative crystallinity and the distribution of spherulite size in core layer with different initial temperatures.

Thus, if the designer wants to improve efficiency, he should cool down the melt with a relative lower initial temperature.

4. Conclusion

We have extended our previous proposed multiscale model and multiscale algorithm to the simulation of 3D polymer crystallization during cooling stage. With the multiscale model and multiscale algorithm, we obtained the temperature distributions and relative crystallinity at various locations in the mold cavity, meanwhile, we also predicted the crystal morphology development and its size as well as its distributions. We have also shown the importance of coupling between the heat transfer with crystallization as well as 3D numerical studies. Results presented in this paper shows that latent heat released by crystallization plays a very important role in macro-micro-coupling which should be considered in order to predict the more accurate crystal morphology; 2D simulation is qualitatively agreement with the 3D simulation (not only for the variables predicted in macroscopic level and microscopic level, but also for the effects of cooling rate and initial temperature), however, in the view of quantitative analysis, results obtained for these two cases have some differences. Therefore, if the computational conditions are allowed, we recommend the 3D simulation.

Future work will concentrate on 3D crystallization simulation during cooling coupled with heat transfer in reinforced system as well as flow-induced crystallization (FIC) simulation in injection molding.

Acknowledgment

The financial supports provided by the Henan Scientific and Technological Research Project (no. 122102210198), Key Scientific and Technological Research Project of Department of Education of Henan Province (no. 12B110006), Youth Scientific Foundation of Henan University of Science and Technology (no. 2012QN015), and the Doctoral Foundation of Henan University of Science and Technology (no. 09001612) are fully acknowledged.