M. Asanishi, T. Takaki and Y. Tomita / Proceedings, AES – ATEMA’ 2007 International Conference,Montreal, Canada, August 06 -10, 2007, pp. 195 - 203

ISBN 0-9780479

International Conference on Advances and Trends in Engineering Materials and their Applications(AES – ATEMA’ 2007)

Polymer Spherulite Growth Simulation
during Crystallization by Phase-Field Method
M. Asanishi1*, T. Takaki2 and Y. Tomita3
1

Graduate School of Science and Technology, Kobe University, Kobe, Japan(Email: asanishi@solid.mech.kobe-u.ac.jp)
2
Faculty of Maritime Sciences, Kobe University, Kobe, Japan(Email: takaki@maritime.kobe-u.ac.jp)
3Faculty of Engineering, Kobe University, Kobe, Japan(Email: tomita@mech.kobe-u.ac.jp)*Corresponding Author

Abstract
The establishment of a coupled numericalmodel that enables us to simulate spherulite
formations and their mechanical behavior
continuously for crystalline polymers is ourultimate goal. In this paper, we focus ourattention on the formation process of spherulite
and the establishment of a spherulite growthmodel. As a numerical model, the phase-field
method is employed here. In the crystallization
process of polymers, viscosity and supercoolingplay an important role in the formation of
microstructures. Here, we introduce theviscosity effect as the degree of molecularalignment at the solid-liquid interface. In thismodel, we employ three order parameters:© AES-Advanced Engineering Solutions (Ottawa, Canada)
All rights are reserved

phase field, which distinguishes between solids
or liquids; crystallographic orientation, whichcharacterizes the growth direction of lamellae;
and temperature. The effects of viscosity andsuper cooling on the formed microstructuralmorphology are investigated by performing aseries of numerical simulations.

Keywords
Spherulite, Solidification, Phase-Field Method,
Polymer, Lamella

1 Introduction
To enable the wide use of semi-crystalline
polymers as structural materials, it isindispensable to establish a numerical modelthat can predict the macroscopic mechanicalM. Asanishi, T. Takaki and Y. Tomita / AES – ATEMA’ 2007 International Conference,Montreal, Canada, August 06 -10, 2007, pp. 195 - 203

behavior taking their microstructures into
consideration. Tomita and Uchida [1], who aremembers of our group, have proposed amultiscale model relating the micro- tomesoscopic
deformation
behaviorsofsemi-crystalline polymers in 2D problems.
Furthermore, they have generalized their 2Dmodel to a 3D model in order to represent the
3D structure of lamellae within a spherulite [2].
Although, as mentioned above, sophisticatedmultiscale models have been developed, theinitial spherulite morphology with a radialarrangement of thin lamellae is the input dataand is assumed when performing thesimulation. The deformation properties ofstructural polymers are intimately related to the
crystallized morphology [3], which changes
from a single crystal to the spherulite structurewith increasing supercooling [4], and thespherulite is classified into different types
depending on the crystallization temperature[5]. Therefore, in order to perform a morecomprehensive investigation considering thedifference in the microstructure, it is crucial toestablish a coupled numerical model of themechanical properties and the morphologicalformation.
Surprisingly few numerical studies have so far
been carried out on spherulite formationconsidering the microstructure, although there
have been many studies on crystallization
kinetics based on the JMAK model [6,7].
Raabe [8] performed a mesoscale simulation of
spherulite growth using a cellular automaton.
However, the results were limited to theevaluation of crystallization kinetics. Kyu et al.
[9, 10] developed phase-field models forcrystalline polymers with simple shapes. Veryrecently, Mattozzi et al. [11] simulated theformation of a spherulite-like structure by a196
© AES-Advanced Engineering Solutions (Ottawa, Canada)
All rights are reserved

lattice Monte Carlo based a random walk.
Although this model includes lengthening,branching and twisting of the lamellae and avisually accurate 3D spherulite structure is
reproduced, the model includes many
assumptions. Granasy et al. [12] successfully
applied a phase-field model for polycrystalline
solidification to the simulation of various
complex morphologies including spherulite.
However, their results were restricted tometallic alloys.
Our ultimate goal is to construct a coupled
numerical model that can simulate themorphology formation and the mechanicalproperties continuously. In this paper, wederive a phase-field model for polymer
spherulites by generalizing the model proposedby Granasy et al. [12].

2 Spherulite
Figure 1 shows the hierarchical structures in apolymer. The spherulite shown in (b) iscomposed of fibrous ribbons growing outward
from the center along the b-axis and its
branching or secondary nucleation. The fibrousribbon (c) is a two-phase material, consistingof a crystalline lamella and an amorphous layer,that grows with twisting and branching. The
lamella has dominant (fold surface) andslower-growing fronts (lateral surface), as
illustrated in (d).
Although using the multiscale model formechanical behavior [2], it is possible to treat
all the structures shown in Fig. 1, in this study,
we aim to treat the structures shown in Figs. 1(a) and (b) because of the restriction to 2D. The
anisotropic interface energy is closely related
to spherulite formation [13]. To simulate thespherulite morphology, the anisotropic growthand branching of lamellae are introduced to theM. Asanishi, T. Takaki and Y. Tomita / AES – ATEMA’ 2007 International Conference,Montreal, Canada, August 06 -10, 2007, pp. 195 - 203

numerical model presented here.

f (I , T )p (I ) f sol (T )
 1  p (I ) f liq (T )  Wq(I ),Figure.1 Hierarchical structures in polymer,(a) crystallizing texture, (b) spherulite,(c) twisting lamella and (d) chain folding
structure.

3 Phase-Field Model
In order to develop the phase-field model ofthe super cooling solidification of purepolyethylene, we develop our model from thefollowing free energy functional [14]:F2ª D (\ )º³«I  f (I , T )  g (I ) s T » dA, (1)¬ 2¼

where I is a phase field that equals zero in
the liquid phase and unity in the solid phase, T
is temperature, T is the normalized
crystallographic orientation with the interval [0,1] and D \ is the gradient coefficientdepending on the angle \ between the interface
normal and the x-axis. f(I7 is the free energy
density in bulk polyethylene, and we employ anormal double well type function:(2)

where W is the barrier height, fsol (T) and fliq (T)are the free energy densities inside the solidand liquid respectively, and we use the freeenergy difference fsol(T)-fliq(T)=L(T-Tm)/Tm as adriving force. p(I and q(I are the energy
interpolating function and the double wellfunction, respectively, and we employp(I  I(10-15II) and q(I  I(I).
In Eq. 1, the third term represents the excess
free energy density due to the presence ofmisorientation at the solid-liquid interface andthe grain boundary. This term enables us tosimulate the branching phenomenon of the
fibrous ribbon or secondary nucleation. Here, s
is a positive constant that determines theamount of excess energy due to themisorientation. g(I is a monotonicallyincreasing function satisfying g(1) = 1 andg(0) = 0 to eliminate the misorientation effectin the liquid phase, and we chooseg(I  I(10-15II).
The interface energy between the fold surfaceof the lamella and the super cooled melt isalmost one order greater than that of the lateralsurface [13]. Such a large difference in(3)interface energies is modelled by adopting thefollowing equations for the gradientcoefficient:

D (\ )
where­D 1 (\ )°®
°D (\ )
¯ 2(for 2Sik+Tm d \ST/kd 2S i+1)k+Tm) (for 2Sik-Tm d \ST/kd 2Sik +Tm),

197© AES-Advanced Engineering Solutions (Ottawa, Canada)
All rights are reserved

M. Asanishi, T. Takaki and Y. Tomita / AES – ATEMA’ 2007 International Conference,Montreal, Canada, August 06 -10, 2007, pp. 195 - 203­ §D ª
2ST ·½º¸ ¾»«1  J cos®k ¨\ 
k ¹ ¿¼
1 J ¬¯ ©

D1 (\ )
D 2 (\ )D1 T m  T
cos \  T .
cos T m

(4)

V

(5)

V

Here, Jis the strength of anisotropy, k is theanisotropic mode, D is a constant and Tm isthe first missing orientation angle. If J k2-1),
D \  D \ for all \
The time evolution equations for ITand T are
obtained as follows:
wI
wtª w § wD wI · w § wD wI ·¸¸  ¨¨ D¸¸M I « ¨¨ D¬ wx © w\ wy ¹ wy © w\ wx ¹º
wf (I , T ) wg (I )s T »   D I 
wI
wI
¼

(6)

2

wT
wt

ª
T º
M T   « g (I ) s
»
T »¼
«¬

2D \
W

2 tanh 1 1  2O ,

W

1
s'T ,
2

(10)

6D \ L(Tm  T )
MI ,
Tm
2W

(11)

3 2



where O is a minimum value of I in theinterface region, such as O 0Althoughfrom Eq. 10, the interface energy changes withthemisorientation

angle

ຒ

T,

the

misorientation dependence on the interface
energy Vis eliminated by setting a small valueof s. Then, from Eq. 7, the constant s relates
only to the rate of reorientation at thesolid-liquid interface.

4 Numerical Simulation(7)wT
L wp (I ) wIN 2T ,               (8)wt
c wI wt
where t is time, Nis the thermal diffusivity andc is the specific heat. MI and MT are themobilities of the phase field I and orientation T,
respectively. The mobility MI determines the
crystallization rate and is derived from thekinetic theory [8]. The mobility MT is related to
the rotational diffusion, which is a function ofviscosity. Both mobilities depend on
temperature, and the ratio MTMI determines
the degree of second nucleation, or branching.
By considering the 1D equilibrium condition at
T = Tm and the steady growth condition, we canobtain the following relations for the interface
thickness G, interface energy V and the constantinterface velocity V:

G

D\

(9)

198
© AES-Advanced Engineering Solutions (Ottawa, Canada)
All rights are reservedFirst, the effects of super cooling and viscosity
on the solidification morphology are
investigated by performing the simulation ofsuper cooling solidification. Since, however,
the growth rate of the polymer is known to bemuch slower than that of thermal diffusion, tocarry out the super cooling solidification
simulation with high accuracy, we need a verylarge computational time. Therefore, using the
hypothetical material constants, we perform aqualitative investigation. Second, in order toperform a more realistic simulation of polymer
spherulite formation, we conduct thesimulation under an isothermal condition.
We use an explicit finite difference scheme tosolve the governing equations. Periodic
boundary conditions are applied to a squarecomputational region. The initial distributionsof the orientation are random in the range [0, 1]
in both the liquid and the nucleus [12].

M. Asanishi, T. Takaki and Y. Tomita / AES – ATEMA’ 2007 International Conference,Montreal, Canada, August 06 -10, 2007, pp. 195 - 203

5 Super-cooling Solidification
5.1 Effects of super-cooling and viscosity
In order to qualitatively investigate the effectsof super cooling and viscosity, i.e., the ratio
MTMI, on the morphology. We perform thesuper cooling solidification simulation using
assumed material constants. Here, we use thefollowing parameters:'x ='y = 10 nm,
G 'x = 40 nm,V= 0.3 J/m2, J 0.03, D =morphology, as shown in Fig. 3 (a), regardlessof the anisotropic interface energy.
It can be concluded from the above numericalresults that the shapes of the morphology andmicrostructures within the spherulite arelargely dependent on the ratio of thecrystallization rate to the rotational one, i.e.,
the ratio MTMI.

1.28010-4 (J/m)0.5, W = 9.888107 J/m3 andMI= 6.327 m3/Js. The computational domain is
set to be 44 Pm = 400 'x400'y (400400grid).
Figure 2 indicates the effect of the driving
force, or super cooling, on the morphology forMTMI = 10. Here, ' is the degree of
normalized super cooling defined by (Tm Ti)/Tm. In the case of small ', the crystal
morphology becomes dendritic in shape, as
shown in Fig. 2 (a), because of the
diffusion-controlled process. In the case oflarge ', it is observed that the crystal is circular,
as can been seen in Figs. 2 (b) and (c).
Furthermore, from Fig. 2, we can see finer
microstructures with increasing driving force,
because the crystallization rate becomes fasterthan the rotational one.
Figure 3 shows the effect of the ratio MTMI on
the morphology for ' = 1.0. It is observed thatthicker and longer crystals grow within
spherulite with increasing MTMI In otherwords, when MTMI becomes smaller than acertain critical point, the growing front has
difficulties in incorporating the same
crystallographic orientation and generates asecond nucleation. The disordered distributionsof the crystallographic orientation along thegrowing front cause the formation of a circular

Figure.2 Effects of driving force onmorphology and microstructure for ' =(a) 0.8 (t = 0.32 Psec), (b) 0.9 (t = 0.19 Psec)
and (c) 1.0 (t = 0.15 Psec).

Figure.3 Effects of mobility ratio on
morphology and microstructure at t = 0.15 Psec
for MTMI = (a) 100, (b) 200 and (c) 300.
5.2 Polymer spherulite
In this section, we discuss the typical spherulite
formation of the polymer, as shown in Fig. 4.
Experimental studies [15] have shown that theearly stage of polymer spherulite growth is theformation of a bundle of lamellae, as shown in
Fig. 4. The growth fronts that induce secondary
nucleation lamellae split and branch. The
branching of lamellae brings about theconversion into ribbon-shaped crystals, such as199© AES-Advanced Engineering Solutions (Ottawa, Canada)
All rights are reserved

M. Asanishi, T. Takaki and Y. Tomita / AES – ATEMA’ 2007 International Conference,Montreal, Canada, August 06 -10, 2007, pp. 195 - 203

those in Fig. 4. The morphology eventually
becomes that of the spherulite by the space
being filled with the branched lamellae.

Figure.4

Spherulite growth sequence.

Figure 5 illustrates the crystal growth in thecase of the high anisotropic interface energy ofJ 0.87he anisotropic mode is employed as
k = 2 to represent the lamellae growth. As
shown in Fig. 5, when a tip of a growing
lamella absorbs the orientational defects into
the growth front, it changes its growing
direction. Furthermore, it can be observed thatas the ratio MTMIdecreases, the branching at
the liquid-solid interface occurs more
frequently and results in the polycrystalline
morphology with fine grains.
On the basis of the surface nucleation model[13] proposed by Hoffman, the growth mode of
the polymer is categorized into three regimes,
i.e., regimes I, II and III. Each regime isdetermined by the ratio of the secondary
nucleation rate to the crystal growth rate at the
chain folding surface. With increasing rate, atransition from regime I to regime III occurs.
Under the condition in regime I, the lamellae
preferentially grow as shown in Fig. 5 (c). In
regime III, a lot of secondary nuclei are
generated and disordered crystal growth isobserved, as shown in Fig. 5 (a). Although it
seems that the condition shown in Fig. 5 (b) isthat of regime II, the calculated morphology isdifferent from that shown in Fig. 4. Next, we
perform the same simulation for MTMI =230, but the computational domain is changed
to 66 Pm (600600 grid). In the early growth

Figure.5 Calculated morphologies with highinterface anisotropy for MTMI = (a) 100,(b) 200 and (c) 300. The left diagram is thephase field, and the right is orientation.

200© AES-Advanced Engineering Solutions (Ottawa, Canada)
All rights are reserved

process, the crystal is thread- shaped, as shownin Fig. 6. Secondary nucleation brings about
splitting and branching lamellae, and finallythe spherulite shown in Fig. 4 is formed. It is
concluded that the formation process of thispolymer spherulite requires both theneedle-shaped lamellae growth and thesecondary nucleation or, branching.

M. Asanishi, T. Takaki and Y. Tomita / AES – ATEMA’ 2007 International Conference,Montreal, Canada, August 06 -10, 2007, pp. 195 - 203

6 Isothermal Solidification
In the previous section, the spherulite
formation simulation of the supercooling melt
is performed using the assumed parameters toqualitatively investigate the effects of thesolidification conditions on the formed
morphology. In this section, in order to enable
a more realistic simulation of the spherulite
formation of the polymer, we perform thesimulation under isothermal conditions inwhich the experimental studies were carriedout [5]. By assuming isothermal conditions, it
is not necessary to solve Eq. 8, and thecomputational effort can be reduced. We use
the following parameters: V= 0.1 J/m2, 'x ='y = 10 nm, G 4'x = 40 nm, D = 7.39010-5(J/m)0.5, W = 3.296107 J/m3, Tm = 418.7 K [13],and MI = 4.087  10-7 m3/sJ. The value ofMIis determined by assuming V = 0.1 Pm/s
[8] and T = 412 K from Eq. 11. Here, thenumerical results for isotropic interface energy
are indicated.
Figure. 7 shows the effect of the mobility ratioon the morphology and microstructure. It isfound that we obtain equivalent results to thoseshown in Fig. 3.
Figure 8 shows a time sequence showing the
growth of the number of spherulites. The
domain is set to be 1212 Pm (12001200

Figure.6 The occurrence of spherulite with
branching lamellae at t = (a) 0.036, (b) 0.99,(c) 0.144, (d) 0.198 and (e) 0.252 Psec.

grid) and MTMI= 60. Twelve nuclei are
generated simultaneously initially (sitesaturated nucleation model). All nuclei growand impinge on the neighbor spherulites.
Finally, the whole space is filled with thespherulite morphologies. Although these results
are for the isotropic interface energy, the201© AES-Advanced Engineering Solutions (Ottawa, Canada)
All rights are reserved

M. Asanishi, T. Takaki and Y. Tomita / AES – ATEMA’ 2007 International Conference,Montreal, Canada, August 06 -10, 2007, pp. 195 - 203

obtained morphology is very similar toexperimentally observed results, e.g., Fig. 1 ofRef. [16].

Figure.7 Spherulite morphology at t = 25 s
for MTMI= (a) 20, (b) 40 and (c) 60

under the isothermal condition was carried out
using realistic material parameters of thepolymer, which provides a spherulite
morphology similar to experimentally observed
results.
Acknowledgments
This research was partially supported by the
Ministry of Education, Culture, Sports, Scienceand Technology Grant-in-Aid for Scientific
Research (B), 18360061, 2007.

References

Figure.8 Spherulite growth sequences t =(a) 5, (b) 15, (c) 25, (d) 35, (e) 45 and (f) 55 s.

7 Conclusions
A phase-field model for polymer spherulite
growth has been established by generalizing
the model proposed by Granasy et al. [12]. By
performing a series of simulations, it isclarified that the driving force of solid-liquidinterface migration and the viscosity, which
control the crystallization and rotation of
crystal orientation, are dominant factors that
determine the morphology and microstructure.
Furthermore, it is confirmed that, in order toreproduce the typical spherulite growth process,a high anisotropy, which creates theneedle-shaped crystals and the secondary
nucleation are required. Finally, a simulation

202© AES-Advanced Engineering Solutions (Ottawa, Canada)
All rights are reserved[1] Tomita, Y. and Uchida, M. (2005)
Computational Characterization of Microto Mesoscopic Deformation Behavior of
Semicrystalline Polymer, Int. J. Mech. Sci.,
47 (4/5), 687 – 700.
[2] Tomita, Y. and Uchida, M. (2005)
Computational
Characterization
andEvaluation of Deformation Behavior of
Spherulite of High Density Polyethylene in
Mesoscale Domain, Computer Modeling in
Engineering & Sciences, 10 (3), 239 - 248.
[3] Way, J. L., Atkinson, J. R. and Nutting, J.,(1974) The effect of spherulite size on
the fracture morphology of polypropylene,
J. Mater. Sci., 9 (2), 293 - 299.
[4] Magill, J. H. and Plazek, D. J. (1967)
Physical
Properties
ofAromatic
Hydrocarbons. II. Solidification Behavior of
1,3,5-Tri-D-Naphthylbenzene*, J. Chem.
Phys., 46 (10), 3757 - 3769.
[5] Norton, D. R. and Keller, A. (1984) The
spherulitic and lamellar morphology ofmelt-crystallized isotactic polypropylene,
Polymer, 26 (5), 704 - 716.
[6] Koscher, E. and Fulchiron, R. (2002)
Influence of shear on polypropylene
crystallization: morphology development

M. Asanishi, T. Takaki and Y. Tomita / AES – ATEMA’ 2007 International Conference,Montreal, Canada, August 06 -10, 2007, pp. 195 - 203

and kinetics, Polymer, 43 (25), 6931 - 6942.
[7] Hieber, C. A. (1995) Correlations for yhequiescent crystallization kinetics of isotactic
polypropyleneandpoly(ethylene
terephthalate, Polymer, 36 (7), 1455 - 1467.
[8] Raabe, D. (2004) Mesoscale simulation of
spherulite
growthduringpolymercrystallization by use of a cellular
automaton, Acta. Mater., 52 (9), 2653 2664.
[9] Kyu, T., Mehta, R. and Chiu, H. –W. (2000)
Spatiotemporal growth of faceted andcurved single crystals, Phys. Rev. E, 61 (4),4161 - 4170.
[10] Metha, R., Keawwattana, W. and Kyu, T.
(2004)
Growth dynamics of isotactic
polypropylene single crystals duringisothermal crystallization from a miscible
polymeric solvent, J. Chem. Phys., 120 (8),4024 - 4031.
[11] Mattozzi, A., Serralunga, P., Hedenqvist,
M. S. and Gedde, U. W. (2006)
Mesoscale modelling of penetrant
diffusion
incomputer-generated
polyethylene–spherulite-likestructures,
Polymer, 47 (15), 5588 - 5595.
[12] Granasy, L., Pusztaiet, T., Tegze, G.,
Warren, J. A. and Douglas, J. F. (2005)
Growth and form of spherulites, Phys. Rev.
E, 72 (1), 011605.
[13] Hoffman, J. D. and Miller, R. J. (1997)
Kinetics of crystallization from the melt
and chain folding in polyethylene
fractions revisited: theory and experimant,
Polymer, 38 (13), 3135 - 3212.
[14] Kobayashi, R., Warren, J. A. and Carter, W.
C. (2000) A continuum model of grainboundaries, Physica D, 140 (1-2), 141 150.

[15] Lei, Y. –G., Chan, C. –M., Wang, Y., Ng,
K. –M., Jiang, Y. and Lin, L. (2003)
Growth process of homogeneously andheterogeneously nucleated spherulites as
observed by atomic force microscopy,
Polymer, 44 (16), 4673 - 4679.
[16] Magill, H. (2001) Spherulites: A personal
perspective*, J. Mater. Sci., 36 (13), 3143
– 3164.

203© AES-Advanced Engineering Solutions (Ottawa, Canada)
All rights are reserved