- Email: [email protected]
Numerical simulation of the polymer crystallization during cooling stage by using level set method
Computational Materials Science 97 (2015) 245–253
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Numerical simulation of the polymer crystallization during cooling stage by using level set method Z.J. Liu, J. Ouyang ⇑, W. Zhou, X.D. Wang Department of Applied Mathematics, Northwestern Polytechnical University, Xi’an 710129, PR China
a r t i c l e
i n f o
Article history: Received 2 May 2014 Received in revised form 20 September 2014 Accepted 20 October 2014 Available online 12 November 2014 Keywords: Simulation Level set Kinetics Morphology Crystallization
a b s t r a c t Aiming to reduce computation complexity in studying crystal interaction by a level set method, we present a new method: the level set method by solving a single signed distance function with the aid of ‘‘colors’’. In this method, the crystal is identified by its assigned color (respectively a number). The growth of a crystal can be captured by semi-Lagrangian method with the aid of coloring. The problem of evolving multiple crystal interfaces is reduced to track one level set variable (signed distance function) and determine the color of a newly solidification node point. The method is used to study the evolution of polymer crystallization under isothermal condition and during cooling stage with a multi-scale model which is proposed in this paper. Good agreement of the numerical results with the theoretical values and those predicted in literature is obtained. The method is easy to be implemented. Compared with the multiple level sets method, it has low computational complexity and low memory requirement. The method is also applicable for any system that displays nucleation and growth. Ó 2014 Elsevier B.V. All rights reserved.
1. Introduction The crystalline structure of semi-crystalline polymers is a major research interest in the field of polymer physics because the physical, mechanical and optical properties of the finished products of polymers depend very strongly on the microstructure [1]. As the demand for complex plastic products rise, it becomes increasingly important to understand and correctly predict morphological formation. Different approaches have been proposed for morphological modeling of polymer crystallization, such as the front-tracking method [2,3], the cellular automata method [4], and the phasefield method [5–7]. In fact, most of the aforementioned methods were presented to simulation of solidification and growth processes for metal systems at first. And they were extended to study the polymer crystallization shortly thereafter. A review of these methods as applied to solidification processes is given in [8]. The level set method is used widely in various applications such as two-phase flow, crack propagation, computer vision and image processing. It has been shown to be a promising mathematical tool for tracking interface with low computation [8]. Numerous studies on crystalline morphology for metal solidification processes have been conducted with level set methods. However, the studies with level set methods for polymer crystallization are still rare.
⇑ Corresponding author. Tel.: +86 29 8849 5234; fax: +86 29 8849 1000. E-mail address: [email protected] (J. Ouyang). http://dx.doi.org/10.1016/j.commatsci.2014.10.038 0927-0256/Ó 2014 Elsevier B.V. All rights reserved.
Moreover, even for metal solidification processes, the efforts for realistic multiple dendrite simulation by using level set methods are not intensified. Ref. [8] indicates that the main reason for this is the huge computation cost in studying crystal interaction with level-set method. The objective of this article is to reduce computation complexity in studying crystal interaction by using level set method. For this purpose, we present a new approach: the level set method by solving a single signed distance function with the aid of ‘‘colors.’’ In this method, the crystal is identified by its assigned color; the growth of a crystal can be captured by semi-Lagrangian method with the aid of coloring; the problem of evolving multiple crystal interfaces is reduced to track one level set variable (signed distance function) and determine the color of a newly crystallization node point. The present paper is organized as follows. Section 2 presents the level set method for the morphological evolution of crystals. In Section 3, the multi-scale model and the coupled computational method are introduced. Section 4 gives the numerical results and discussions. Finally the conclusions are drawn in Section 5.
2. The level set method for polymer crystallization The level set method was first proposed by Osher and Sethian [9] for front propagation with curvature dependent speed. Since then, it has been extended to numerous applications with moving
246
Z.J. Liu et al. / Computational Materials Science 97 (2015) 245–253
(a) Two crystals are identified by its node color at time t n-1
(b) Uncolored nodes in the crystalline phase are captured by level set method at time t n
(c) Color the nodes by semi-Lagrangian method at time tn Fig. 1. Schematic representation of the level set method for polymer crystallization.
Fig. 2. Semi-Lagrangian transport procedure.
interfaces in fluid mechanics, combustion, computer animation, image processing, among others. As a novel approach of handling the interface evolution, it implied in a manner to express the closed planar curves, three-dimensional closed surface or, more generally, interface in Rn of co-dimension one, so as to avoid the evolution process of tracing the interface. Thus problem of interface evolution is converted to a pure searching for the numerical answers of partial differential equation. In fact, in the level set method the interface is represented as the zero contour of a level set function which satisfies an evolution equation. The physics which governs the motion of the interface motion is included naturally in this model [10]. A level set function /(x, t) is defined as:
8 > < þdðx; tÞ x 2 X /ðx; tÞ ¼ 0 x 2 @ X ¼ CðtÞ > : dðx; tÞ x 2 X
ð2:1Þ
where C(t) is an interface in Rn of co-dimension one, X is an open region bounded by C(t), and d(x, t) is the smallest distance between a given point in the domain and the interface:
dðx; tÞ ¼ min ðjx xC j; jy yC jÞ xC ;yC 2C
ð2:2Þ
The interface is equal to the zero level set of /, i.e.
CðtÞ ¼ ðx 2 X : /ðx; tÞ ¼ 0Þ
ð2:3Þ
Z.J. Liu et al. / Computational Materials Science 97 (2015) 245–253
The normal direction of the interface is calculated as follows:
n¼
r/ jr/j
ð2:4Þ
The evolution equation of /, corresponding to the motion of the face, is given by the convection equation:
/t þ ujr/j ¼ 0
ð2:5Þ
After an update according to the level set Eq. (2.5), / does not in general remain a signed distance function. It is thus necessary for re-initialization where the following equation is iterated until reaching steady-state:
/0 ð1 jr/jÞ /t ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi /20 þ e2
ð2:6Þ
where /0 is the initial level set value to be re-initialized. After / reaches steady-state, |r/| = 1, i.e. / is a signed distance. The parameter e in Eq. (2.6) takes some small value and is needed for the formulation to remain well-posed as / ? 0. When tracking only a crystal, it is convenient to use a level set to indicate the position of the interface. However it is not a very easy task to determine the boundaries when multiple crystals make contact and form boundaries. Multiple level sets method is used to solve this problem [11,12]. But solving multiple level set equations is only realistic for a small number of crystals. As the number of crystals increases, the multiple level sets method becomes increasingly inefficient. Similarly to [8] we use a single signed distance to implicitly represent the interface of crystals. In order to identify different crystals, we assign each crystal a ‘‘color’’ (respectively a number). As demonstrated in Fig. 1, the dotted lines and the solid line are the interfaces of the two crystals, which are captured by level set method at time tn1 and tn (tn = tn1 + Dt). The nodes inside the two crystals are colored blue (the left) and red (the right) respectively (see Fig. 1(a)). The nodes inside the solid line and outside the dotted line are in the crystalline phase at time tn, they are captured by level set method and the colors of them need to be determined (see Fig. 1(b)). Noting Semi-Lagrangian transport keep certain material properties along the trajectory of
247
material point, it naturally emerges as an attractive method to determine colors of the nodes. The trajectory that passes through the node at time tn can be followed in time to tn1 to obtain the old location. The color of the node is the same of the old location. If the old location is colored node, there is no problem; if not, it is obvious we cannot determine the color of the old location by interpolation. Fortunately, there is a simple alternative that we can set the color of the old location to the color of the nearest colored node. The impingement boundary of the two crystals is shown by the colored nodes (see Fig. 1(c)). Fig. 2 illustrates the SemiLagrangian transport procedure. The nodes inside two crystals are colored red and blue respectively, and the nodes whose colors need to be determined are gray (Fig. 2(a)). Each gray nodal point is moved backwards in time (Fig. 2(b)) by using the classical Runge– Kutta midpoint method [13]:
x ¼ xn uðxn ; t n ÞDt=2
ð2:7Þ
xn1 ¼ xn uðx ; t n ÞDt
ð2:8Þ
At the old location, the color is determined by the nearest colored node (Fig. 2(c)), leading to the color at the nodal point xn (Fig. 2(d)). The validity of the level set method for studying crystal interaction is shown in Section 4. 3. Modeling for polymer crystallization During the cooling stage, the process of crystallization is coupled with the heat transfer. Any realistic simulation of polymer crystallization should not only account for the nucleation, growth and impingement of spherulites at the microscopic level, but also for the heat diffusion at the macroscopic level. Referring to [14], we use multi-scale model to simulate the process of polymer crystallization. 3.1. Thermal field in the macroscopic scale The energy equation plays a role in macro–micro coupling by a latent heat release term. The corresponding simplified energy equation is [14]:
qcp
@T @a ¼ r ðkp rTÞ þ qDH @t @t
ð3:1Þ
where q is the material density, cp is the thermal capacity, t is the time, kp is the thermal conductivity, DH is the total enthalpy released during crystallization, and a is the relative crystallinity. The last source term is the contribution of the latent heat released by the crystallization. Crystal Temperature Mesoscale
3.2. Crystal morphology evolution in the microscopic scale Crystal microstructure evolution involves the processes of nucleation, growth and impingement of spherulites. We simulate the nucleation of spherulites using the following relation of nucleation density [3,15]:
n h io2=3 NðTÞ ¼ 1:458 N0 exp u T 0m T
Coarse grid
Fine grid
Fig. 3. Schematic representation of multi-scale algorithm.
ð3:2Þ
where T 0m is the equilibrium melting temperature, N0 and u are the material dependent parameters. A stochastic method is used in this study to place the nucleation sites in the case of isothermal crystallization. A new nucleus appears randomly and occupies a node, i.e. a node is chosen randomly, if the node is colored, delete the nucleus; if not, color the node by a new color and update the signed distance field with the following operation:
248
Z.J. Liu et al. / Computational Materials Science 97 (2015) 245–253
Fig. 4. Multi-scale algorithm.
(a) t =0s
(c) t =80s
(b) t =32s
(d) Final morphology
Fig. 5. Examples of polymer morphology obtained from the Level set method simulation at Tc = 390 K.
Z.J. Liu et al. / Computational Materials Science 97 (2015) 245–253
249
the interface. In this paper, u denotes the normal speed of growth (the growth rate), and it is first computed at the nodal points near the interface by (3.4), then at other points by a algorithm to extend the interface velocity away from the interface as described in [8]. It should be noted that a node is marked as being near the interface if at least one of its neighboring nodes has a different sign of /.
3.3. Macro–micro-algorithm
Fig. 6. Relative crystallinity evolution under different temperature in isothermal condition.
uðyÞ ¼ maxðu0 ðyÞ; kx yk R0 Þ; 8y 2 X
ð3:3Þ
where x is the nucleation sites, /0 is the signed distance before the potential nucleation site is nucleated at x, R0 is the size of the initial crystal seed at location x and y is the location of a node. It is worthwhile to note that the algorithm for placement of new nuclei can be extended to the non-isothermal case easily, more details about it is given in Section 3.3. We simulate the growth and impingement of spherulites based on the method developed in Section 2. According to [16], the growth rate of the spherulites is given by:
GðTÞ ¼ G0 exp
" # U Kg exp m RðT T 1 Þ T T0 T
Our multi-scale model for polymer crystallization incorporates two distinct length scales to simulate the microstructure: a coarse grid for the heat diffusion, a fine grid for tracking the solid fraction evolution. Referring to [14], we build up the algorithm as the following: The domain is first divided by a coarse grid. Then, each coarse grid is subdivided into a number of cells to form the fine grid. This is schematically shown in Fig. 3. On the coarse grid, the FVM is used to solve the energy equation. On the fine grid the Level set method is implemented to capture the growth front of a population of spherulites. In our model, we assign a color (respectively a number) to each crystal and then assign this color (number) to the spatial locations (nodes) covered by the crystal. Therefore, the relative crystallinity of each coarse grid can be calculated by
ð3:4Þ
where G0 and Kg are constants, U⁄ is the activation energy of motion, R is the gas constant, T1 is a temperature typically 30 K below the glass transition. In the level set method, the interface velocity u should be defined on the whole domain or on a narrow band near
a¼
number of nodes that have been colored total number of nodes
ð3:5Þ
It is through (3.1) and (3.5) that the microstructural level of nucleation and growth and the macroscopic level of heat diffusion are coupled. It is worth noting that we suppose the temperature is the same in a control volume of the coarse grid in non-isothermal case. i.e., we do not consider the spatial inconsistency of the temperature in a control volume. Therefore, on one hand, we can place new nuclei in a coarse grid in non-isothermal case by using the algorithm in the isothermal case. In this way, if temperature changes, no modification is needed for the algorithm, only the number of the nuclei to be placed need a corresponding change according to the temperature in a control volume; on the other hand, It is easier to present the algorithm which is used for extending the interface velocity away from the freezing interface to every other nodal point in a coarse grid. Because the temperature of every nodal point is the same in a control volume, we give the same speed to each nodal point in a coarse grid according to (3.4). The macro–micro-algorithm used for our simulation is shown in Fig. 4. In the algorithm, the index j denotes the variable at time tj, and i denotes the variable in the ith control volume.
4. Results and discussions 4.1. Problem formulation
Fig. 7. Comparison at constant cooling rate of relative crystallinity between Kolmogoroff–Avrami model and level set method.
Numerical experiment is performed on a rectangle of size 4 mm 4 mm. We cool the sample starting from a uniform temperature of T0 and then cooling the boundary with constant speed c. The material we consider here is isotactic polypropylene (iPP) and parameters used in the simulation are [17] N0 = 17.4 106/ m3, u = 0.155, v0 = 2.1 104 m/s, U⁄/Rg = 755 k, Kg = 534858 K2, T 0m ¼ 467 K, Tg = 266 K, q = 900 kg/m3, cp = 2.14 103 J/(kg k), kp = 0.193 W/(m K), DHX1 = 107 103 J kg. Unless otherwise stated, the other parameters are set to be T0 = 470 K, C = 10 K/min. The number of coarse grid is 6 6 and the number of fine grid in a coarse grid is 150 150.
250
Z.J. Liu et al. / Computational Materials Science 97 (2015) 245–253
(a)
(c)
(b)
(e)
(d)
(f)
Fig. 8. Evolution of crystal morphology in a control volume. (a) t = 138 s. (b) t = 178 s. (c) t = 202 s. (d) t = 218 s. (e) t = 232 s. (f) Final morphology.
4.2. Validity of the level set method for crystallization In order to validate the Level set method for crystallization, we compare the numerical results with the theoretical values in the isothermal and non-isothermal conditions respectively. 4.2.1. Isothermal case In the case of isothermal crystallization, the number of nuclei N and the growth rate G for a given temperature Tc are calculated from (3.2) and (3.4), The calculated number of nuclei appears instantaneously and randomly within a square lattice with 600 600 grids at the onset of crystallization (t = t0 = 0), The growth front of the spherulites is captured by the level set method for crystallization. Fig. 4 illustrates the morphological development during the isothermal crystallization of iPP at Tc = 390 K obtained from the level set method for crystallization, the size of the rectangle is 200 lm 200 lm. Fig. 5(a) shows iPP morphology at t = 0 s when heterogeneous nucleation occurs. Up until t = 32 s, each spherulite grows individually without impingement (see Fig. 5(b)). As time increases, several impingement boundaries occur (see Fig. 5(c)) and finally the final morphology is observed (see Fig. 5(d)). Avrami [18,19] proposed an isothermal model which has been used universally to describe polymer crystallization kinetics. The change in crystallinity with time can be readily expressed as:
aðtÞ ¼ 1 expðka tna Þ
ð4:1Þ
where a is the relative crystallinity at time t; na, the Avrami index (crystal geometry information); and ka, the isothermal crystallization rate constant containing the nucleation and growth rates. For our case, na = 2, ka = pN(Tc)[G(Tc)]2 [20]. Thus, Eq. (4.1) can be rewritten as:
aðtÞ ¼ 1 expfpNðT c Þ½GðT c Þ2 t2 g
ð4:2Þ
Fig. 6 shows the effect of temperature on the overall crystallization kinetics in isothermal conditions. Symbols show the predicted results from the simulation, and lines show the analytical values according to Avrami model. The simulated results match well with the analytical values. This clearly confirms the validity of the level set method for crystallization. 4.2.2. Non-isothermal case In the case of non-isothermal polymer crystallization, we here present the results obtained in a control volume (a coarse grid). It is worth noting that the temperature used here is T = T0 ct which is independent of evolution of relative crystallinity. In the non-isothermal polymer crystallization, Kolmogoroff– Avrami model [15,21] is one of the important models in description of polymer crystallization kinetics. In this model, Kolmogoroff model describes crystallinity evolution accounting of the number of nuclei per unit volume and spherulitic growth rate. The corresponding equations of the Kolmogoroff–Avrami model are:
a ¼ 1 expðaf Þ af ¼ C m
Z 0
t
dNðsÞ ds
ð4:3Þ Z
t
m GðuÞdu ds
ð4:4Þ
s
where a is the relative crystallinity, af is the fictive volume fraction, m is a dimensionality exponent, Cm is a shape factor, G is the linear growth rate and N is the nucleation density. In our case, Cm = p, m = 2. The numerical results of crystallization kinetics compare to the values calculated with the Kolmogoroff–Avrami model which is shown in Fig. 7. As we can see that the numerical results show a good agreement with the theoretical values. It should be noted that the accuracy of the level set method for crystallization is dependent on the grid size: the finer grid size, the more accurate of the computational results. Moreover, the less time step will also lead to better results. However, more CPU time is needed for these
Z.J. Liu et al. / Computational Materials Science 97 (2015) 245–253
(a)
(b)
(c)
(e)
251
(d)
(f)
Fig. 9. Evolution of crystal morphology in the whole computation zone. (a) t = 120 s. (b) t = 169 s. (c) t = 190 s. (d) t = 216 s. (e) t = 246 s. (f) Final morphology.
methods. In this simulation, the time step is set as Dt = 0.1 s. Fig. 8 shows the morphological development under constant cooling rate. It is obvious that most of the impingement boundaries are curves while they are lines in the isothermal conditions. The structure and characteristics of crystals obtained in this study accorded with the literature [3] Thus, the level set method for crystallization presented here is valid and reliable. 4.3. The results of multi-scale algorithm Fig. 9 shows the evolution of crystal morphology in the whole computation zone. In the figure, the white region is polymer melt,
the other colored region is the crystals. It shows clearly the effect of the boundary cooling on the development of the morphology. Most crystals nucleate and grow fast close to the boundary. The generated morphology is typical for an experiment with boundary cooling, where there are many fine grains close to the boundary and relative small amounts of larger ones inside. Fig. 10 shows the evolution of relative crystallinity, it is clear that the relative crystallinity of the control volumes close to the boundary reach 1 first. Fig. 11 shows the evolution of temperature; it is obvious that the relative crystallinity change in accordance with the temperature of control volume. Evolution of temperature at polymer core for different cooling rates is shown in Fig. 12, we can find that a
252
Z.J. Liu et al. / Computational Materials Science 97 (2015) 245–253
(a) t=15s
(b) t =197s
(c) t =214s
(d) t=224s
(e) t =239s
(f) t =282s
Fig. 10. Evolution of relative crystallinity.
(a) t =15s
(b) t =197s
(c) t =214s
(d) t =224s
(e) t =239s
(f) t =282s
Fig. 11. Evolution of temperature.
Z.J. Liu et al. / Computational Materials Science 97 (2015) 245–253
253
5. Conclusions
Fig. 12. Temperature evolution at polymer core for different cooling rates.
We have presented an efficient level set method to simulate the solidification microstructures of crystalline and semicrystalline polymer. The procedure involves using ‘‘color’’ to identify different crystals and determining the color of the nodes by signed distance function and semi-Lagrangian method. We also have built up a multi-scale model for the crystallization during the cooling stage: a coarse grid for the heat diffusion and a finer grid for microstructure evolution. To solve this multi-scale model, we have presented an efficient multi-scale algorithm. We also have shown the effects of cooling rates. In conclusion, we have simulated polymer crystallization under isothermal and non-isothermal conditions by level set method. It is easy to be implemented. And it is applicable for any system that displays nucleation and growth because the motion of the interface in the level set method is determined by a velocity field regardless of the shape of the interface. Compared with the multiple level sets method, the method has low computational complexity and low memory requirement. What is more, the predicted data by our simulation scheme for polymer crystallization show fair agreement with the theoretical values, it illustrates that our scheme can not only present the detail information about the microstructure evolution of polymer crystallization but also provide the reliable crystallization kinetics data. It may be a useful tool for the optimal designs of the molding process. Acknowledgment All the authors would like to acknowledge support by the National Basic Research Program of China (2012CB025903). References
Fig. 13. Relative crystallinity in the whole computation zone for different cooling rates.
quasi-isothermal plateau is developed in different cooling rates. And it is clear that the larger cooling rate leads to an earlier show of the temperature plateau. Fig. 13 shows relative crystallinity in the whole computation zone for different cooling rates, it is obvious the lower cooling rate leads to a longer time for crystallization process.
[1] S. Anantawaraskul, S. Ketdee, P. Supaphol, J. Appl. Polym. Sci. 111 (2009) 2260– 2268. [2] S. Swaminarayan, C. Charbon, Polym. Eng. Sci. 38 (1998) 634–643. [3] C. Charbon, S. Swaminarayan, Polym. Eng. Sci. 38 (1998) 644–656. [4] D. Raabe, A. Godara, Modell. Simul. Mater. Sci. Eng. 13 (2005) 733–751. [5] M. Asanishi, T. Takaki; Y. Tomita, Proceedings of AES-ATEMA’2007 International Conference, Montreal, Canada, 2007, pp. 195–203. [6] H. Xu, C.T. Bellehumeur, J. Appl. Polym. Sci. 102 (2006) 5903–5917. [7] H. Xu, C.T. Bellehumeur, J. Appl. Polym. Sci. 107 (2008) 236–245. [8] L.J. Tan, Cornell University (PhD Thesis), Ithaca, 2007. [9] S. Osher, J.A. Sethian, J. Comput. Phys. 79 (1988) 12–49. [10] X. Zhang, J.S. Chen, S. Osher, Inter. Multi. Mech. 1 (2008) 178–191. [11] M.O. Bloomfield, D.F. Richards, T.S. Cale, Philos. Mag. 83 (2003) 3549–3568. [12] G. Russo, P. Smereka, SIAM J. Sci. Comput. 21 (2000) 2073–2095. [13] D.R. Durran, Numerical Methods for Fluid Dynamics: With Applications to Geophysics, Springer, 2010. [14] C. Ruan, J. Ouyang, S. Liu, Int. J. Heat. Mass. Transfer. 55 (2012) 1911–1921. [15] R. Pantani, I. Coccorullo, V. Speranza, G. Titomanlio, Prog. Polym. Sci. 30 (2005) 1185–1222. [16] J.D. Hoffman, R.L. Miller, Polymer 38 (1997) 3151–3212. [17] B. Yang, X.R. Fu, W. Yang, L. Huang, M.B. Yang, J.M. Feng, Polym. Eng. Sci. 48 (2008) 1707–1717. [18] M. Avrami, J. Chem. Phys. 7 (1939) 1103–1112. [19] M. Avrami, J. Chem. Phys. 8 (1940) 212–224. [20] J.M. Schultz, Polymer Crystallization: The Development of Crystalline Order in Thermoplastic Polymers, American Chemical Society, Washington, DC, 2001. [21] R. Zheng, R.I. Tanner, X.J. Fan, Injection Molding: Integration of Theory and Modeling Methods, Springer, Berlin, 2011.
Contents lists available at ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Numerical simulation of the polymer crystallization during cooling stage by using level set method Z.J. Liu, J. Ouyang ⇑, W. Zhou, X.D. Wang Department of Applied Mathematics, Northwestern Polytechnical University, Xi’an 710129, PR China
a r t i c l e
i n f o
Article history: Received 2 May 2014 Received in revised form 20 September 2014 Accepted 20 October 2014 Available online 12 November 2014 Keywords: Simulation Level set Kinetics Morphology Crystallization
a b s t r a c t Aiming to reduce computation complexity in studying crystal interaction by a level set method, we present a new method: the level set method by solving a single signed distance function with the aid of ‘‘colors’’. In this method, the crystal is identified by its assigned color (respectively a number). The growth of a crystal can be captured by semi-Lagrangian method with the aid of coloring. The problem of evolving multiple crystal interfaces is reduced to track one level set variable (signed distance function) and determine the color of a newly solidification node point. The method is used to study the evolution of polymer crystallization under isothermal condition and during cooling stage with a multi-scale model which is proposed in this paper. Good agreement of the numerical results with the theoretical values and those predicted in literature is obtained. The method is easy to be implemented. Compared with the multiple level sets method, it has low computational complexity and low memory requirement. The method is also applicable for any system that displays nucleation and growth. Ó 2014 Elsevier B.V. All rights reserved.
1. Introduction The crystalline structure of semi-crystalline polymers is a major research interest in the field of polymer physics because the physical, mechanical and optical properties of the finished products of polymers depend very strongly on the microstructure [1]. As the demand for complex plastic products rise, it becomes increasingly important to understand and correctly predict morphological formation. Different approaches have been proposed for morphological modeling of polymer crystallization, such as the front-tracking method [2,3], the cellular automata method [4], and the phasefield method [5–7]. In fact, most of the aforementioned methods were presented to simulation of solidification and growth processes for metal systems at first. And they were extended to study the polymer crystallization shortly thereafter. A review of these methods as applied to solidification processes is given in [8]. The level set method is used widely in various applications such as two-phase flow, crack propagation, computer vision and image processing. It has been shown to be a promising mathematical tool for tracking interface with low computation [8]. Numerous studies on crystalline morphology for metal solidification processes have been conducted with level set methods. However, the studies with level set methods for polymer crystallization are still rare.
⇑ Corresponding author. Tel.: +86 29 8849 5234; fax: +86 29 8849 1000. E-mail address: [email protected] (J. Ouyang). http://dx.doi.org/10.1016/j.commatsci.2014.10.038 0927-0256/Ó 2014 Elsevier B.V. All rights reserved.
Moreover, even for metal solidification processes, the efforts for realistic multiple dendrite simulation by using level set methods are not intensified. Ref. [8] indicates that the main reason for this is the huge computation cost in studying crystal interaction with level-set method. The objective of this article is to reduce computation complexity in studying crystal interaction by using level set method. For this purpose, we present a new approach: the level set method by solving a single signed distance function with the aid of ‘‘colors.’’ In this method, the crystal is identified by its assigned color; the growth of a crystal can be captured by semi-Lagrangian method with the aid of coloring; the problem of evolving multiple crystal interfaces is reduced to track one level set variable (signed distance function) and determine the color of a newly crystallization node point. The present paper is organized as follows. Section 2 presents the level set method for the morphological evolution of crystals. In Section 3, the multi-scale model and the coupled computational method are introduced. Section 4 gives the numerical results and discussions. Finally the conclusions are drawn in Section 5.
2. The level set method for polymer crystallization The level set method was first proposed by Osher and Sethian [9] for front propagation with curvature dependent speed. Since then, it has been extended to numerous applications with moving
246
Z.J. Liu et al. / Computational Materials Science 97 (2015) 245–253
(a) Two crystals are identified by its node color at time t n-1
(b) Uncolored nodes in the crystalline phase are captured by level set method at time t n
(c) Color the nodes by semi-Lagrangian method at time tn Fig. 1. Schematic representation of the level set method for polymer crystallization.
Fig. 2. Semi-Lagrangian transport procedure.
interfaces in fluid mechanics, combustion, computer animation, image processing, among others. As a novel approach of handling the interface evolution, it implied in a manner to express the closed planar curves, three-dimensional closed surface or, more generally, interface in Rn of co-dimension one, so as to avoid the evolution process of tracing the interface. Thus problem of interface evolution is converted to a pure searching for the numerical answers of partial differential equation. In fact, in the level set method the interface is represented as the zero contour of a level set function which satisfies an evolution equation. The physics which governs the motion of the interface motion is included naturally in this model [10]. A level set function /(x, t) is defined as:
8 > < þdðx; tÞ x 2 X /ðx; tÞ ¼ 0 x 2 @ X ¼ CðtÞ > : dðx; tÞ x 2 X
ð2:1Þ
where C(t) is an interface in Rn of co-dimension one, X is an open region bounded by C(t), and d(x, t) is the smallest distance between a given point in the domain and the interface:
dðx; tÞ ¼ min ðjx xC j; jy yC jÞ xC ;yC 2C
ð2:2Þ
The interface is equal to the zero level set of /, i.e.
CðtÞ ¼ ðx 2 X : /ðx; tÞ ¼ 0Þ
ð2:3Þ
Z.J. Liu et al. / Computational Materials Science 97 (2015) 245–253
The normal direction of the interface is calculated as follows:
n¼
r/ jr/j
ð2:4Þ
The evolution equation of /, corresponding to the motion of the face, is given by the convection equation:
/t þ ujr/j ¼ 0
ð2:5Þ
After an update according to the level set Eq. (2.5), / does not in general remain a signed distance function. It is thus necessary for re-initialization where the following equation is iterated until reaching steady-state:
/0 ð1 jr/jÞ /t ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi /20 þ e2
ð2:6Þ
where /0 is the initial level set value to be re-initialized. After / reaches steady-state, |r/| = 1, i.e. / is a signed distance. The parameter e in Eq. (2.6) takes some small value and is needed for the formulation to remain well-posed as / ? 0. When tracking only a crystal, it is convenient to use a level set to indicate the position of the interface. However it is not a very easy task to determine the boundaries when multiple crystals make contact and form boundaries. Multiple level sets method is used to solve this problem [11,12]. But solving multiple level set equations is only realistic for a small number of crystals. As the number of crystals increases, the multiple level sets method becomes increasingly inefficient. Similarly to [8] we use a single signed distance to implicitly represent the interface of crystals. In order to identify different crystals, we assign each crystal a ‘‘color’’ (respectively a number). As demonstrated in Fig. 1, the dotted lines and the solid line are the interfaces of the two crystals, which are captured by level set method at time tn1 and tn (tn = tn1 + Dt). The nodes inside the two crystals are colored blue (the left) and red (the right) respectively (see Fig. 1(a)). The nodes inside the solid line and outside the dotted line are in the crystalline phase at time tn, they are captured by level set method and the colors of them need to be determined (see Fig. 1(b)). Noting Semi-Lagrangian transport keep certain material properties along the trajectory of
247
material point, it naturally emerges as an attractive method to determine colors of the nodes. The trajectory that passes through the node at time tn can be followed in time to tn1 to obtain the old location. The color of the node is the same of the old location. If the old location is colored node, there is no problem; if not, it is obvious we cannot determine the color of the old location by interpolation. Fortunately, there is a simple alternative that we can set the color of the old location to the color of the nearest colored node. The impingement boundary of the two crystals is shown by the colored nodes (see Fig. 1(c)). Fig. 2 illustrates the SemiLagrangian transport procedure. The nodes inside two crystals are colored red and blue respectively, and the nodes whose colors need to be determined are gray (Fig. 2(a)). Each gray nodal point is moved backwards in time (Fig. 2(b)) by using the classical Runge– Kutta midpoint method [13]:
x ¼ xn uðxn ; t n ÞDt=2
ð2:7Þ
xn1 ¼ xn uðx ; t n ÞDt
ð2:8Þ
At the old location, the color is determined by the nearest colored node (Fig. 2(c)), leading to the color at the nodal point xn (Fig. 2(d)). The validity of the level set method for studying crystal interaction is shown in Section 4. 3. Modeling for polymer crystallization During the cooling stage, the process of crystallization is coupled with the heat transfer. Any realistic simulation of polymer crystallization should not only account for the nucleation, growth and impingement of spherulites at the microscopic level, but also for the heat diffusion at the macroscopic level. Referring to [14], we use multi-scale model to simulate the process of polymer crystallization. 3.1. Thermal field in the macroscopic scale The energy equation plays a role in macro–micro coupling by a latent heat release term. The corresponding simplified energy equation is [14]:
qcp
@T @a ¼ r ðkp rTÞ þ qDH @t @t
ð3:1Þ
where q is the material density, cp is the thermal capacity, t is the time, kp is the thermal conductivity, DH is the total enthalpy released during crystallization, and a is the relative crystallinity. The last source term is the contribution of the latent heat released by the crystallization. Crystal Temperature Mesoscale
3.2. Crystal morphology evolution in the microscopic scale Crystal microstructure evolution involves the processes of nucleation, growth and impingement of spherulites. We simulate the nucleation of spherulites using the following relation of nucleation density [3,15]:
n h io2=3 NðTÞ ¼ 1:458 N0 exp u T 0m T
Coarse grid
Fine grid
Fig. 3. Schematic representation of multi-scale algorithm.
ð3:2Þ
where T 0m is the equilibrium melting temperature, N0 and u are the material dependent parameters. A stochastic method is used in this study to place the nucleation sites in the case of isothermal crystallization. A new nucleus appears randomly and occupies a node, i.e. a node is chosen randomly, if the node is colored, delete the nucleus; if not, color the node by a new color and update the signed distance field with the following operation:
248
Z.J. Liu et al. / Computational Materials Science 97 (2015) 245–253
Fig. 4. Multi-scale algorithm.
(a) t =0s
(c) t =80s
(b) t =32s
(d) Final morphology
Fig. 5. Examples of polymer morphology obtained from the Level set method simulation at Tc = 390 K.
Z.J. Liu et al. / Computational Materials Science 97 (2015) 245–253
249
the interface. In this paper, u denotes the normal speed of growth (the growth rate), and it is first computed at the nodal points near the interface by (3.4), then at other points by a algorithm to extend the interface velocity away from the interface as described in [8]. It should be noted that a node is marked as being near the interface if at least one of its neighboring nodes has a different sign of /.
3.3. Macro–micro-algorithm
Fig. 6. Relative crystallinity evolution under different temperature in isothermal condition.
uðyÞ ¼ maxðu0 ðyÞ; kx yk R0 Þ; 8y 2 X
ð3:3Þ
where x is the nucleation sites, /0 is the signed distance before the potential nucleation site is nucleated at x, R0 is the size of the initial crystal seed at location x and y is the location of a node. It is worthwhile to note that the algorithm for placement of new nuclei can be extended to the non-isothermal case easily, more details about it is given in Section 3.3. We simulate the growth and impingement of spherulites based on the method developed in Section 2. According to [16], the growth rate of the spherulites is given by:
GðTÞ ¼ G0 exp
" # U Kg exp m RðT T 1 Þ T T0 T
Our multi-scale model for polymer crystallization incorporates two distinct length scales to simulate the microstructure: a coarse grid for the heat diffusion, a fine grid for tracking the solid fraction evolution. Referring to [14], we build up the algorithm as the following: The domain is first divided by a coarse grid. Then, each coarse grid is subdivided into a number of cells to form the fine grid. This is schematically shown in Fig. 3. On the coarse grid, the FVM is used to solve the energy equation. On the fine grid the Level set method is implemented to capture the growth front of a population of spherulites. In our model, we assign a color (respectively a number) to each crystal and then assign this color (number) to the spatial locations (nodes) covered by the crystal. Therefore, the relative crystallinity of each coarse grid can be calculated by
ð3:4Þ
where G0 and Kg are constants, U⁄ is the activation energy of motion, R is the gas constant, T1 is a temperature typically 30 K below the glass transition. In the level set method, the interface velocity u should be defined on the whole domain or on a narrow band near
a¼
number of nodes that have been colored total number of nodes
ð3:5Þ
It is through (3.1) and (3.5) that the microstructural level of nucleation and growth and the macroscopic level of heat diffusion are coupled. It is worth noting that we suppose the temperature is the same in a control volume of the coarse grid in non-isothermal case. i.e., we do not consider the spatial inconsistency of the temperature in a control volume. Therefore, on one hand, we can place new nuclei in a coarse grid in non-isothermal case by using the algorithm in the isothermal case. In this way, if temperature changes, no modification is needed for the algorithm, only the number of the nuclei to be placed need a corresponding change according to the temperature in a control volume; on the other hand, It is easier to present the algorithm which is used for extending the interface velocity away from the freezing interface to every other nodal point in a coarse grid. Because the temperature of every nodal point is the same in a control volume, we give the same speed to each nodal point in a coarse grid according to (3.4). The macro–micro-algorithm used for our simulation is shown in Fig. 4. In the algorithm, the index j denotes the variable at time tj, and i denotes the variable in the ith control volume.
4. Results and discussions 4.1. Problem formulation
Fig. 7. Comparison at constant cooling rate of relative crystallinity between Kolmogoroff–Avrami model and level set method.
Numerical experiment is performed on a rectangle of size 4 mm 4 mm. We cool the sample starting from a uniform temperature of T0 and then cooling the boundary with constant speed c. The material we consider here is isotactic polypropylene (iPP) and parameters used in the simulation are [17] N0 = 17.4 106/ m3, u = 0.155, v0 = 2.1 104 m/s, U⁄/Rg = 755 k, Kg = 534858 K2, T 0m ¼ 467 K, Tg = 266 K, q = 900 kg/m3, cp = 2.14 103 J/(kg k), kp = 0.193 W/(m K), DHX1 = 107 103 J kg. Unless otherwise stated, the other parameters are set to be T0 = 470 K, C = 10 K/min. The number of coarse grid is 6 6 and the number of fine grid in a coarse grid is 150 150.
250
Z.J. Liu et al. / Computational Materials Science 97 (2015) 245–253
(a)
(c)
(b)
(e)
(d)
(f)
Fig. 8. Evolution of crystal morphology in a control volume. (a) t = 138 s. (b) t = 178 s. (c) t = 202 s. (d) t = 218 s. (e) t = 232 s. (f) Final morphology.
4.2. Validity of the level set method for crystallization In order to validate the Level set method for crystallization, we compare the numerical results with the theoretical values in the isothermal and non-isothermal conditions respectively. 4.2.1. Isothermal case In the case of isothermal crystallization, the number of nuclei N and the growth rate G for a given temperature Tc are calculated from (3.2) and (3.4), The calculated number of nuclei appears instantaneously and randomly within a square lattice with 600 600 grids at the onset of crystallization (t = t0 = 0), The growth front of the spherulites is captured by the level set method for crystallization. Fig. 4 illustrates the morphological development during the isothermal crystallization of iPP at Tc = 390 K obtained from the level set method for crystallization, the size of the rectangle is 200 lm 200 lm. Fig. 5(a) shows iPP morphology at t = 0 s when heterogeneous nucleation occurs. Up until t = 32 s, each spherulite grows individually without impingement (see Fig. 5(b)). As time increases, several impingement boundaries occur (see Fig. 5(c)) and finally the final morphology is observed (see Fig. 5(d)). Avrami [18,19] proposed an isothermal model which has been used universally to describe polymer crystallization kinetics. The change in crystallinity with time can be readily expressed as:
aðtÞ ¼ 1 expðka tna Þ
ð4:1Þ
where a is the relative crystallinity at time t; na, the Avrami index (crystal geometry information); and ka, the isothermal crystallization rate constant containing the nucleation and growth rates. For our case, na = 2, ka = pN(Tc)[G(Tc)]2 [20]. Thus, Eq. (4.1) can be rewritten as:
aðtÞ ¼ 1 expfpNðT c Þ½GðT c Þ2 t2 g
ð4:2Þ
Fig. 6 shows the effect of temperature on the overall crystallization kinetics in isothermal conditions. Symbols show the predicted results from the simulation, and lines show the analytical values according to Avrami model. The simulated results match well with the analytical values. This clearly confirms the validity of the level set method for crystallization. 4.2.2. Non-isothermal case In the case of non-isothermal polymer crystallization, we here present the results obtained in a control volume (a coarse grid). It is worth noting that the temperature used here is T = T0 ct which is independent of evolution of relative crystallinity. In the non-isothermal polymer crystallization, Kolmogoroff– Avrami model [15,21] is one of the important models in description of polymer crystallization kinetics. In this model, Kolmogoroff model describes crystallinity evolution accounting of the number of nuclei per unit volume and spherulitic growth rate. The corresponding equations of the Kolmogoroff–Avrami model are:
a ¼ 1 expðaf Þ af ¼ C m
Z 0
t
dNðsÞ ds
ð4:3Þ Z
t
m GðuÞdu ds
ð4:4Þ
s
where a is the relative crystallinity, af is the fictive volume fraction, m is a dimensionality exponent, Cm is a shape factor, G is the linear growth rate and N is the nucleation density. In our case, Cm = p, m = 2. The numerical results of crystallization kinetics compare to the values calculated with the Kolmogoroff–Avrami model which is shown in Fig. 7. As we can see that the numerical results show a good agreement with the theoretical values. It should be noted that the accuracy of the level set method for crystallization is dependent on the grid size: the finer grid size, the more accurate of the computational results. Moreover, the less time step will also lead to better results. However, more CPU time is needed for these
Z.J. Liu et al. / Computational Materials Science 97 (2015) 245–253
(a)
(b)
(c)
(e)
251
(d)
(f)
Fig. 9. Evolution of crystal morphology in the whole computation zone. (a) t = 120 s. (b) t = 169 s. (c) t = 190 s. (d) t = 216 s. (e) t = 246 s. (f) Final morphology.
methods. In this simulation, the time step is set as Dt = 0.1 s. Fig. 8 shows the morphological development under constant cooling rate. It is obvious that most of the impingement boundaries are curves while they are lines in the isothermal conditions. The structure and characteristics of crystals obtained in this study accorded with the literature [3] Thus, the level set method for crystallization presented here is valid and reliable. 4.3. The results of multi-scale algorithm Fig. 9 shows the evolution of crystal morphology in the whole computation zone. In the figure, the white region is polymer melt,
the other colored region is the crystals. It shows clearly the effect of the boundary cooling on the development of the morphology. Most crystals nucleate and grow fast close to the boundary. The generated morphology is typical for an experiment with boundary cooling, where there are many fine grains close to the boundary and relative small amounts of larger ones inside. Fig. 10 shows the evolution of relative crystallinity, it is clear that the relative crystallinity of the control volumes close to the boundary reach 1 first. Fig. 11 shows the evolution of temperature; it is obvious that the relative crystallinity change in accordance with the temperature of control volume. Evolution of temperature at polymer core for different cooling rates is shown in Fig. 12, we can find that a
252
Z.J. Liu et al. / Computational Materials Science 97 (2015) 245–253
(a) t=15s
(b) t =197s
(c) t =214s
(d) t=224s
(e) t =239s
(f) t =282s
Fig. 10. Evolution of relative crystallinity.
(a) t =15s
(b) t =197s
(c) t =214s
(d) t =224s
(e) t =239s
(f) t =282s
Fig. 11. Evolution of temperature.
Z.J. Liu et al. / Computational Materials Science 97 (2015) 245–253
253
5. Conclusions
Fig. 12. Temperature evolution at polymer core for different cooling rates.
We have presented an efficient level set method to simulate the solidification microstructures of crystalline and semicrystalline polymer. The procedure involves using ‘‘color’’ to identify different crystals and determining the color of the nodes by signed distance function and semi-Lagrangian method. We also have built up a multi-scale model for the crystallization during the cooling stage: a coarse grid for the heat diffusion and a finer grid for microstructure evolution. To solve this multi-scale model, we have presented an efficient multi-scale algorithm. We also have shown the effects of cooling rates. In conclusion, we have simulated polymer crystallization under isothermal and non-isothermal conditions by level set method. It is easy to be implemented. And it is applicable for any system that displays nucleation and growth because the motion of the interface in the level set method is determined by a velocity field regardless of the shape of the interface. Compared with the multiple level sets method, the method has low computational complexity and low memory requirement. What is more, the predicted data by our simulation scheme for polymer crystallization show fair agreement with the theoretical values, it illustrates that our scheme can not only present the detail information about the microstructure evolution of polymer crystallization but also provide the reliable crystallization kinetics data. It may be a useful tool for the optimal designs of the molding process. Acknowledgment All the authors would like to acknowledge support by the National Basic Research Program of China (2012CB025903). References
Fig. 13. Relative crystallinity in the whole computation zone for different cooling rates.
quasi-isothermal plateau is developed in different cooling rates. And it is clear that the larger cooling rate leads to an earlier show of the temperature plateau. Fig. 13 shows relative crystallinity in the whole computation zone for different cooling rates, it is obvious the lower cooling rate leads to a longer time for crystallization process.
[1] S. Anantawaraskul, S. Ketdee, P. Supaphol, J. Appl. Polym. Sci. 111 (2009) 2260– 2268. [2] S. Swaminarayan, C. Charbon, Polym. Eng. Sci. 38 (1998) 634–643. [3] C. Charbon, S. Swaminarayan, Polym. Eng. Sci. 38 (1998) 644–656. [4] D. Raabe, A. Godara, Modell. Simul. Mater. Sci. Eng. 13 (2005) 733–751. [5] M. Asanishi, T. Takaki; Y. Tomita, Proceedings of AES-ATEMA’2007 International Conference, Montreal, Canada, 2007, pp. 195–203. [6] H. Xu, C.T. Bellehumeur, J. Appl. Polym. Sci. 102 (2006) 5903–5917. [7] H. Xu, C.T. Bellehumeur, J. Appl. Polym. Sci. 107 (2008) 236–245. [8] L.J. Tan, Cornell University (PhD Thesis), Ithaca, 2007. [9] S. Osher, J.A. Sethian, J. Comput. Phys. 79 (1988) 12–49. [10] X. Zhang, J.S. Chen, S. Osher, Inter. Multi. Mech. 1 (2008) 178–191. [11] M.O. Bloomfield, D.F. Richards, T.S. Cale, Philos. Mag. 83 (2003) 3549–3568. [12] G. Russo, P. Smereka, SIAM J. Sci. Comput. 21 (2000) 2073–2095. [13] D.R. Durran, Numerical Methods for Fluid Dynamics: With Applications to Geophysics, Springer, 2010. [14] C. Ruan, J. Ouyang, S. Liu, Int. J. Heat. Mass. Transfer. 55 (2012) 1911–1921. [15] R. Pantani, I. Coccorullo, V. Speranza, G. Titomanlio, Prog. Polym. Sci. 30 (2005) 1185–1222. [16] J.D. Hoffman, R.L. Miller, Polymer 38 (1997) 3151–3212. [17] B. Yang, X.R. Fu, W. Yang, L. Huang, M.B. Yang, J.M. Feng, Polym. Eng. Sci. 48 (2008) 1707–1717. [18] M. Avrami, J. Chem. Phys. 7 (1939) 1103–1112. [19] M. Avrami, J. Chem. Phys. 8 (1940) 212–224. [20] J.M. Schultz, Polymer Crystallization: The Development of Crystalline Order in Thermoplastic Polymers, American Chemical Society, Washington, DC, 2001. [21] R. Zheng, R.I. Tanner, X.J. Fan, Injection Molding: Integration of Theory and Modeling Methods, Springer, Berlin, 2011.