166 Macromol. Theory Simul.

9, 166–175 (2000)

Full Paper: This work employs the relaxed Stefan model slowly with time, and finally, the crystallinity of the inter-
and Nakamura crystallization kinetics to describe the non- nal part exceeds the crystallinity on the surface. In the sec-
isothermal crystallization process of polymeric materials ond case (for instance, in water), crystallinity is relatively
by finite element discretization method (FEM) simulation, low, and there is a serious gradient of crystallinity. The
giving the evolution of crystallinity distribution on 2 D crystallinity on the surface reaches a very low equilibrium
space. Numerical results show that the final crystallinity value in a short time and changes little afterwards.
and its distribution are mainly dependent on the cooling Although the crystallinity of the inside part can be
rate. Crystallinity decreases with increasing cooling rate, improved by changing the shape of the polymeric article,
but the influence is negligible as long as the cooling rate the crystallinity on the surface essentially remains con-
is below a critical value (ca. 30 8C N min–1 for poly(ethyl- stant, which leads to a significant gradient. Geometrical
ene terephthalate) (PET)). If the cooling rate is higher shape and dimension of the article are also important to
than this critical one, crystallinity drops sharply. It is also the crystallinity and its distribution, and the ratio of sur-
concluded that the crystallization behavior of polymeric face area to volume can be used as a rough index to esti-
samples in a mild cooling medium is quite different from mate the shape/dimension influence on crystallinity.
that in a strong cooling medium. In the first case (for Except the coefficient of thermal conductivity, physical
example, in silicon oil), crystallinity of the article is rela- parameters of the polymeric material and kinetic para-
tively high and its distribution is fairly uniform. During meters of crystallization show only weak effects com-
the initial short period, the crystallinity on the surface is pared to cooling conditions.
higher than that on the inside. Crystallinity increases

FEM simulation of nonisothermal crystallization, 1

Crystallinity distribution on 2D space
Deyue Yan* 1, Hong Jiang1, Huxi Li2
1
College of Chemistry and Chemical Technology, Shanghai Jiao Tong University, 800 Dongchuan Road, Shanghai 200240,
P. R. China
2
Laboratoire des Sciences du Génie Chimique (LSGC, CNRS UPR 6811), Institute National Polytechnique de Lorraine
(INPL), 1, rue Grandville, B. P. 451, 54001 Nancy Cedex, France
(Received: February 8, 1999; revised: August 9, 1999)

Introduction Ozawa8), and Berger9) models, have been developed over

On processing a crystalline polymer, crystallinity and the years in order to apply them to nonisothermal situa-
morphology are subject variations in different parts of the tions. It has been pointed out that none of these models
product as a result of the non-uniform temperature distri- currently available is totally satisfactory in predicting
bution created by rapid cooling/heating on the surface nonisothermal crystallization based on isothermal data.
combined with low thermal conductivity of most poly- So far, the Nakamura model has proved to be the most
meric materials. Microstructure distribution plays a sig- suitable one. Chan and Isayev10) carried out a comprehen-
nificant role in determining the mechanical properties of sive investigation on this model. They concluded that it is
the product. The distribution of crystallinity or micro- possible to obtain a good description of nonisothermal
structure is dependent on thermal-mechanical histories, crystallization kinetics on using the Nakamura equation
crystallization kinetics and other processing parameters. solely based on isothermal differential scanning chroma-
In general, it can be described by overall crystallization tography (DSC) measurements. Flow-induced crystalliza-
kinetics and thermal-mechanical histories, if only the tion is a more complicated process, which has been
crystallinity is concerned, which will be fully discussed advanced in recent years11–13).
in this paper. During polymer processing, crystallization Numerous experimental and simulation works have
usually occurs under nonisothermal conditions. Various been conducted in order to describe the kinetics of crys-
modifications of the Avrami isothermal crystallization tallization, but few of them aim to characterize the distri-
kinetics1–3), for examples, Malkin4, 5), Nakamura6, 7), bution of microstructure. A computer simulation has been

Macromol. Theory Simul. 9, No. 3 i WILEY-VCH Verlag GmbH, D-69451 Weinheim 2000 1022-1344/2000/0303–0166$17.50+.50/0
FEM simulation of nonisothermal crystallization, 1 ... 167

developed to predict the morphology evolution during where s is the time step length. The boundary and initial
isothermal crystallization14–17). For non-isothermal crys- conditions are
tallization, Chan and Isayev10) reported experimental
results of the crystallinity distribution throughout the qn N m ˆ hT n ‡ gN ns† on CN
thickness direction of a quenched slab. Mazzullo et al.
investigated the phenomena responsible for skin-core T 0 ˆ T0 ; w0 ˆ 0; q0 ˆ 0 7†
morphology by FEM simulation18). This work mainly
focuses on the microstructure (including crystallinity dis- Here gN (ns) is a function of time and m the unitary
tribution, crystalline structure and flow-induced orienta- exterior normal vector on the natural boundary, and CN
tion) of articles manufactured from semi-crystalline poly- the natural boundary. T0 is the initial melt temperature, w0
mers in various processing operations. Moreover, this the initial crystallinity, and q0 the initial heat flux. Using
paper puts emphasis on the pattern of dynamic evolution the standard Galerkin method, we can rewrite Eq. (4) as
of the crystallinity distribution for quiescent crystalliza- follows:
tion under different conditions. Z
‰qCp T n ÿ T nÿ1 † ‡ qL wn ÿ wnÿ1 † ‡ sC N qn ŠNdX ˆ 0 8†
X
Mathematical treatment
where N is one of the shape functions and X is the
Finite element discretization of the problem: It is difficult domain under consideration. Substituting equality
to find an analytical solution for nonisothermal crystalli-
zation on multi-dimensional spaces (2D or 3D). There-
sN†C N qn ˆ C N sNqn † ÿ C sN† N qn 9†
fore, finite element method (FEM) and finite differential
method (FDM) were employed to numerically investigate and natural boundary condition qn N m = hTn + gN (ns) into
the nonisothermal crystallization behavior. The energy Eq. (8), we obtain
conservation equation can be described by
qCp T n ; N† ‡ sa kT n ; N† ‡ sphT n ; NP ˆ qCp T nÿ1 ; N†
qT qw
qCp ‡ qL ‡ CNq ˆ 0 1†
qt qt ‡ qL wnÿ1 ÿ wn †; N† ÿ spgN ns†; NP 10†
where q, Cp and L are density, heat capacity, and heat of with
fusion, respectively. For the sake of simplicity they are
Z
assumed to be constant, although in reality, they depend
u; m† ˆ u N mdX
on crystallinity, temperature, pressure and other factors. X
Related discussion was given elsewhere19). Here, w is the
Z !
crystallinity and q is the thermal flux vector, which is Xd
qu qm
described by the constitutive equation: a u; m† ˆ dX
X iˆ1
qx i qxi
q ˆ ÿkCT 2† Z
pu; mP ˆ u N mdC
where k is the thermal conductivity and T is the tempera- CN
ture. The crystallization kinetics was generalized by a dif-
ferential equation: where d is the dimension of the article studied (d = 2 for
two-dimensional space).
qw Eq. (10) is an integral formula and can not be solved
ˆ f T; w† 3†
qt directly. In order to obtain temperature and crystallinity
distributions, domain X is divided into small sub-regions
The exact expression of the right side depends on the
(meshes or elements), for example, triangular meshes.
crystallization kinetics actually employed.
Temperature and crystallinity distributions inside the ele-
Eq. (1) – (3) are coupled with one another. We
ments are uniquely determined by their values at the
employed the semi-explicit FDM algorithm to realize dis-
nodes and the shape functions. In this way, the continuous
cretization of the continuous time variable:
integral formula of Eq. (10) is transformed into a set of
linear equations, which can be solved by Gauss elimina-
T n ÿ T nÿ1 wn ÿ wnÿ1
qCp ‡ qL ‡ C N qn ˆ 0 4† tion or other P
suitable methods.
s s
Let, T n ˆ jˆ1 Tjn Nj , where M is the number of nodes,
M

qn ˆ ÿkCT n 5† Nj is the shape function of the jth node and Tj is the tem-
wn ÿ wnÿ1 perature of the jth node. Substitution of Tn into Eq. (10)
ˆ f T nÿ1 ; wnÿ1 † 6†
s leads to
168 D. Yan, H. Jiang, H. Li

X
M where E is the number of elements. With T supposed to
‰ qCp Nj ; Ni † ‡ sa kNj ; Ni † ‡ sphNj ; Ni PŠTjn ˆ be linear and w, q, k, Cp, L being constants inside the ele-
jˆ1
ments, we have
X
M 2 3
qCp Nj ; Ni †Tjnÿ1 ‡ qL wnÿ1 ÿ wn †; Ni † ÿ spgN ns†; Ni P 0D e† 0q e† Cpe† 2 1 1
jˆ1 K 1e† ˆ 41 2 15 17a†
24
1 1 2
i ˆ 1; 2; 3; :::; M 11†
0D e† 0 k e† T
or in the matrix form: K 2e† ˆ B B 17b†
2

‰K 1† ‡ sK 2† ‡ sK 3† ŠT n ˆ K 1† T nÿ1 ‡ f 1†
ÿ sf 2†
12† h
K 3e† ˆ
6
where
2 3
2 dij lij ‡ dim lim † dij lij dim lim
Z E Z
X X
E
N4 dij lij 2 dij lij ‡ djm ljm † djm ljm 5 17c†
Kij1† ˆ qCp Ni Nj dX ˆ qCp Ni Nj dX ˆ Kij1e† 13a†
X Xe
dim lim djm ljm 2 dim lim ‡ djm ljm †
eˆ1 eˆ1

where
X
E
Kij
2†
ˆ a kNi ; Nj † ˆ K ij
2e†
13b†
eˆ1 xi yi 1
D e† ˆ det xj yj 1 i; j; m ˆ node indices of element e
X
E xm ym 1
Kij3† ˆ phNi ; Nj P ˆ Kij3e† 13c†
eˆ1  
1 yj ÿ ym ym ÿ yi yi ÿ yj
Bˆ
and D e† xm ÿ xj xi ÿ xm xj ÿ xi
n
X
E 1 nodes i, j 7 CN;e
fi d† ˆ fi de† d ˆ 1; 2 14† dij ˆ
0 otherwise
eˆ1

q
Superscript e represents elements, and E is the total
lij ˆ xi ÿ xj † ‡ yi ÿ yj †
2 2

number of elements.
Eq. (12) is the general formula of the discrete FEM and
equations. Its dimension and therefore computation time
depends on the number of nodes. Some comments about 2 3
1
0D e† 0 q e† L e†
Eq. (12) – (14) are given below: (a) K(d) (d = 1, 2, 3) are ˆ wnÿ1 ÿ wn † 4 1 5
1e† e†
f 18a†
M 6 M matrices and f(d) (d = 1, 2) are M 6 1 vectors. In 6
1
our simulation, M varies from 1 200 to 4 800. In other
words, Eq. (12) is a set that is composed of 1 200 – 4 800 2 3
linear equations. Special techniques are required to dij lij gN;i ns† ‡ dim lim gN;i ns†
1 4
reduce the storage requirement for matrices K and to f 2e†
ˆ dij lij gN;j ns† ‡ djm ljm gN;i ns† 5 18a†
2
increase computation speed. For this sake, we employed dim lim gN;m ns† ‡ djm ljm gN;m ns†
sparse matricesP to store K. (b)P In serial addition opera-
E de† E de† where gN,k (t) is the value of gN of the kth node at time t.
tions, such as K
eˆ1 ij and eˆ1 fi , the local node
It is easy to calculate the coefficients of Eq. (12) by
indices must be mapped to be the global node indices in
substituting Eq. (17) and (18) into Eq. (13) and (14). If
advance. The exact expressions for coefficients K and f
physical properties are independent of temperature and
are dependent on the discretization of the computational
crystallinity, K(d) is also independent of temperature and
domain, or in other words, the mesh type. In this paper,
crystallinity, and can be computed. f (d) is the function of
we employ triangular mesh. The two-dimensional domain
temperature and crystallinity and thus, must be updated at
X was decomposed into triangular meshes {Xe}, and
each step. During programming, f (1) should be divided
boundary C was decomposed into segments {CN,e}, or
into two parts. One of them is related to the crystallinity
of the (n-1)th step, which must be computed before updat-
X ˆ [Eeˆ1 Xe 15†
ing the crystallinity, while the other one is related to the
crystallinity of the nth step and must be computed after
CN ˆ [Eeˆ1 CN;e 16† updating the crystallinity.
FEM simulation of nonisothermal crystallization, 1 ... 169

In the substitution operation of Eq. (17) and (18) into Cp, heat conductivity k and latent heat L, are listed as fol-
Eq. (13) and (14), we have to transform the local node lows:
indices into the global node indices. A triangular mesh
composed of three nodes is taken into account below. Let U* = 6 284 J N mol–1 Tg = 73 8C
its global node indices be i, j, and m, and its local node
indices be 1, 2, and 3, respectively, then K11de† and K12de† Tm = 280 8C Tv = Tg – 30 8C
d† d†
hould be added to KP ii and Kij , respectively, in the serial
E de†
addition operation, eˆ1 Kij . This rule is also true for (1/t1/2)0 = 3.58 6 104 s–1 Kg = 3.66 6 105 K2
other components of K(de), and it is also applicable to the
q = 1.1567 g N cm–3 Cp = 2.259 6 107 erg N g–1 K–1
addition operation of f (de). The local node indices can be
1, 2, or 3 for triangular elements. The global node indices k = 1.9 6 104 erg N s–1 N cm–1 N K–1 L = 0.26 6 109erg N g–1
range from 1 to M, where M is the number of total nodes.
Each node has a unique global node index. In fact, the In this paper, we use PET as a model material to ana-
order of the global indices in FEM computation is very lyze the crystallinity distribution under nonisothermal
important to reduce computation time and storage conditions. The material parameters of PET listed above
requirement. There are several methods available to opti- will be employed unless particularly specified. The max-
mize node numbering. imum crystallinity of PET is 0.33.
The discretization of the heat flux constitutive equation
(Eq. (5)) is very simple. For isotropic materials it reads
Results and discussion
 
k e† yj ÿ ym †Tin ‡ ym ÿ yi †Tjn ‡ yi ÿ yj †Tmn The influence of the cooling rate on crystallinity: The
qn ˆ ÿ e† 19†
D xm ÿ xj †Tin ‡ xi ÿ xm †Tjn ‡ xj ÿ xi †Tmn typical dependence of crystallization rate on temperature,
K (T), is shown in Fig. 1. The crystallization rate is sensi-
Nonisothermal crystallization kinetics: Based on iso-
tive to temperature. No crystallization occurs beyond a
thermal kinetic conditions and the assumption that the
certain temperature range, and there is a maximum crys-
number of activated nuclei keeps constant, Nakamura et
tallization rate at a temperature near the middle of the
al.6, 7) developed an overall crystallization kinetics from
range. Crystallinity is determined by the “effective crys-
the Avrami theory. The differential form of the Nakamura
tallization period” during nonisothermal crystallization.
equation is:
The so-called “effective crystallization period” is defined
by a temperature interval, in which crystallization occurs
dh and beyond which no crystallization can be observed. For
ˆ nK T† 1 ÿ h†‰ÿln 1 ÿ h†Š
nÿ1†=n
20†
dt example, the effective crystallization period is very short,
where n is the Avrami index obtained from isothermal and crystallinity is extremely low if the melt is quickly
crystallization experiments, and K(T) is the nonisother- quenched in liquid nitrogen, which is a widely used pro-
mal crystallization rate constant at temperature T, which cedure for freezing-in the microstructure of polymers.
is related to the Avrami isothermal crystallization rate The non-uniform distribution of crystallinity originates
constant by from the difference in cooling rate at different parts of the
manufactured articles. In most examples, different parts
 
1
K T† ˆ ‰ln 2†Š
1=n
21†
t1=2
Here, t1/2 is the half time, and it can be expressed by:
 
1
ˆ
t1=2
     
1 U  =R Kg
exp ÿ exp ÿ 22†
t1=2 0
T ÿ Tv T N DT N f

where R is the universal gas constant, DT = Tm – T is the

supercooling, f = 2T/(T + Tm) is a correction factor
accounting for the reduction in the latent heat of fusion as
the temperature drops, and (1/t1/2)0, U*, and Tv are mate-
rial constants. Kinetic parameters and other material para- Fig. 1. Crystallization rate constant K(T) of PET as a function
meters of PET, such as density q, specific heat capacities of temperature
170 D. Yan, H. Jiang, H. Li

(a)

Fig. 2. The influence of cooling rates on the ultimate crystalli-

nity of PET. The maximum crystallinity (cooling rate
approaches zero) is assumed to be unit
(b)
of an article manufactured under nonisothermal condi-
tions have different thermal histories, and therefore dif-
ferent cooling rates, and different effective crystallization
periods, respectively. In general, the surface of the article
is cooled much more quickly than the inside part during
processing, which results in low crystallinity on the sur-
face and probable deterioration of the surface quality of
articles. Fig. 2 shows ultimate crystallinities of PET sam-
ples uniformly cooled at cooling rates ranging from
1 8C N min–1 to 200 8C N min–1. It can be used as a guide to
estimate the influence of thermal history on the crystalli-
nity. We can expect that the faster the cooling rate, the
shorter the effective crystallization period and the lower
the ultimate crystallinity. For example, the relative crys-
tallinity of the article decreases from 100% to 28% with
increasing the cooling rate from 20 8C N min–1 to Fig. 3. Temperature distribution (a) and crystallinity distribu-
100 8C N min–1. Fig. 2 reveals that crystallinity is almost tion (b) of an infinitely extended PET slab quenched in water at
independent of the cooling rate at low cooling rates 92 8C for 200 s simulated by setting kx = 0 and ky = 1.9 6 104
(a30 8C N min–1 for PET). So it is possible to produce arti- erg N s–1 N cm–1 N K–1
cles with fairly uniform crystallinity distribution by keep-
ing the cooling rate below the critical value. This critical Isayev10) reported experimental results of the crystallinity
cooling rate is dependent on the crystallization kinetics of distribution in an infinitely extended slab with finite
the polymeric material used. For PET, the critical cooling thickness quenched on both sides (i. e., the length is much
rate is ca. 30 8C N min–1. Crystallinity decreases dramati- greater than the width of the rectangular cross-section, Lx
cally if the cooling rate is higher than this critical value. S Ly). There are two possible approaches to simulate the
Reliability of simulation: Computer simulation is a experiment. At first, we can directly simulate temperature
powerful tool to investigate complicated situations which and crystallinity distributions of a rectangular cross-sec-
cannot be examined by experimental approaches due to tion satisfying: Lx S Ly. The difficulty of this method is
economical and technical difficulties. For example, it is that it requires a large number of elements to obtain reli-
very time consuming and costly to measure the crystalli- able precision. In fact this difficulty can be avoided if we
nity at every point of an article in order to obtain the dis- rewrite Eq. (19) into
tribution of crystallinity, but it can be easily carried out
by computer simulation. However, simulation means 1
qn ˆ ÿ
modeling and simplification of real problems. Consider- D e†
ing the possible error introduced by numerical methods   
or program defects (bugs), it is necessary to verify the k 0 yj ÿ ym †Tin ‡ ym ÿ yi †Tjn ‡ yi ÿ yj †Tmn
N x 23†
realibility of a simulation by experimental data. Chan and 0 ky xm ÿ xj †Tin ‡ xi ÿ xm †Tjn ‡ xj ÿ xi †Tmn
FEM simulation of nonisothermal crystallization, 1 ... 171

Fig. 4. Crystallinity distribution of a 3 mm thickness PET sam- Fig. 5. Crystallinity distribution of a 6 mm thickness PET sam-
ple quenched in silicon oil at temperatures of 24, 81 and 120 8C. ple quenched in water at temperatures of 71 and 92 8C. Symbols
Symbols denote experimental data and lines represent model denote experimental data and lines represent model predictions;
predictions; h = 58 J N s–1 N m–2 N K–1, kx = ky = 1.9 6 104 h = 420 J N s–1 N m–2 N K–1, kx = ky = 1.9 6 104 erg N s–1 N cm–1 N K–1
erg N s–1 N cm–1 N K–1

and set kx = 0. This is equivalent to the experimental con-

dition: Lx S Ly. Fig. 3 shows temperature distribution (a)
and crystallinity distribution (b) simulated by setting kx =
0. There is no gradient of temperature or crystallinity
along the X-direction. These are exactly the conditions
on which the experiments were implemented.
A comparison of simulation results with experimental
data of a slab quenched in silicon oil at 24, 81 and
120 8C, respectively, is shown in Fig. 4. The predicted
results fit the measured data very well at all of the three
different cooling temperatures. Both calculated and meas-
ured data of crystallinity show little variation along the
thickness direction even at low cooling temperature
(24 8C). If the thickness of the article increases to 6 mm,
and the sample is quenched in water (a much stronger
cooling medium than silicon oil), a significant gradient of
Fig. 6. Crystallinity distribution after cooling in silicon oil for
crystallinity is observed from both simulation and experi- 40 s at 81 8C; h = 58 J N s–1 N m–2 N K–1, kx = ky = 1.9 6 104
ment as shown in Fig. 5. The heat transfer coefficient erg N s–1 N cm–1 N K–1
employed was 420 J N s–1 N m–2 N K–1, which is consider-
ably higher than the correspondent value of silicon oil.
This leads to quick cooling on the surface and a signifi- is dependent on the shape of the article, material proper-
cant crystallinity gradient. Simulation results were further ties, and cooling environment.
confirmed by the similarity between predicted crystalli- Dealing with a sample which is much longer than the
nity distributions and measured ones. In the simulation, dimension of its cross-section, we can only investigate
only heat transfer coefficient, h, between cooling medium the crystallinity distribution in this cross-section. Fig. 6 –
and sample is the variable, which equals 58 and 420 9 show the evolution of crystallinity distribution in the
J N s–1 N m–2 N K–1 in oil and water, respectively. 3 mm 6 12 mm rectangular cross-section of a sample
Evolution of crystallinity distributions by cooling in a cooled in silicon oil. After cooling for 40 s (Fig. 6), the
mild medium: In Fig. 2, we have demonstrated the rela- crystallinity near the surface reaches a certain degree,
tionship between crystallinity and cooling rate for a uni- while it is still negligible on the inside. Apparently, the
formly cooled sample. In a real processing process, this surface is quickly cooled below the melting point (Tm),
relationship is much more complicated. Here, tempera- and begins to crystallize. Nevertheless, the temperature
tures are non-uniform at different parts inside an article, of the center part is still very high due to the low heat
depending on the interaction of conduction, heat genera- conductivity of polymeric materials and therefore, crys-
tion, and heat exchange with surroundings, which in turn tallinity is very low. For example, the crystallinity on the
172 D. Yan, H. Jiang, H. Li

Fig. 7. Crystallinity distribution after cooling in silicon oil for Fig. 9. Crystallinity distribution after cooling for 200 s in sili-
80 s at 81 8C; h = 58 J N s–1 N m–2 N K–1, kx = ky = 1.9 6 104 con oil at 81 8C; h = 58 J N s–1 N m–2 N K–1, kx = ky = 1.9 6 104
erg N s–1 N cm–1 N K–1 erg N s–1 N cm–1 N K–1

tion, which is only a little different from the distribution

given in Fig. 8. After 200 s, crystallinity does not change
any further, and the temperature in the whole article has
reached ambient temperature. It is interesting that the
final crystallinity distribution along the thickness direc-
tion is of fairly uniform open contours. Concisely, rapid
cooling of the boundary caused that the final crystallinity
on the surface remains at a lower value than the one in
the center part. Silicon oil is a fairly mild cooling med-
ium, and the final crystallinity reaches as high as 50%
which is corresponding to an average cooling rate of
approximatly 60 8C N min–1 (Fig. 2).
Evolution of crystallinity distributions by cooling in a
strong medium: The objective of this part is to investigate
the evolving pattern of crystallinity distribution under
strong cooling conditions and to understand the relation-
Fig. 8. Crystallinity distribution after cooling for 100 s in sili- ship between crystallinity distribution and cooling envi-
con oil at 81 8C; h = 58 J N s–1 N m–2 N K–1, kx = ky = 1.9 6 104
erg N s–1 N cm–1 N K–1 ronments.
The evolution of crystallinity distribution for a sample
with 3 mm 6 12 mm rectangular cross-section which is
surface has reached about 10%, while it is less than 2% at cooled in water is shown in Fig. 10 – 12. Compared to
the center. Crystallinity reaches ca. 30% after cooling for cooling in silicon oil, crystallinity is much lower, and
80 s in silicon oil (Fig. 7). The distribution pattern is quite crystallization is much faster. For example, if crystalliza-
different from that of 40 s. At first, crystallinity increases tion is terminated after 60 s, the maximum crystallinity is
from the surface to the center and reaches a maximum ca. 4% only. The average crystallinity is ca. 3%, indicat-
value. After that it decreases again and reaches a mini- ing an average cooling rate of A 200 8C N min–1 according
mum at the center. This is similar to the classic crystalli- to Fig. 2. The second feature is the large gradient of crys-
zation front, moving from surface to center. Fig. 8 and 9 tallinity distribution. On the boundary, crystallinity is
display crystallinity distributions after cooling for 100 much lower than in the center, with the latter one being
and 200 s, respectively. After 100 s (Fig. 8), the shape of relatively low itself. If a sample is quenched in water, the
the crystallinity distribution is inverted: the center part temperature on the boundary drops rapidly because of the
shows higher crystallinity than the surface, and the crys- strong cooling ability of water, which leads to low crys-
tallinity decreases gradually from the center to the sur- tallinity on the surface. Because of the low thermal con-
face. Fig. 9 indicates the ultimate crystallinity distribu- ductivity of polymeric materials, the internal section
FEM simulation of nonisothermal crystallization, 1 ... 173

Fig. 10. Crystallinity distribution after cooling for 18 s in Fig. 12. Crystallinity distribution after cooling for 60 s in
water at 81 8C; h = 420 J N s–1 N m–2 N K–1, kx = ky = 1.9 6 104 water at 81 8C; h = 420 J N s–1 N m–2 N K–1, kx = ky = 1.9 6 104
erg N s–1 N cm–1 N K–1 erg N s–1 N cm–1 N K–1

Influence of the sample shape on crystallinity and its

distributions: It has been pointed out in the previous sec-
tion that the crystallinity of the 3 mm 6 12 mm rectangu-
lar cross-section is quite different from that of the 6 mm
thickness slab. The maximum crystallinity of the former
is only ca. 4% if quenched in water at 81 8C, while the
latter reaches as high as 30% and 50% if quenched in
water at 71 8C and 92 8C. The significant difference in
crystallinity originates from different cooling histories,
which in turn arise from the geometrical difference. In
order to analyze the effect of geometrical shapes of arti-
cles on the crystallinity and its distribution, we define

Cooling Area
kˆ 24†
Volume
where “cooling area” is the area of the surface that is in
Fig. 11. Crystallinity distribution after cooling for 30 s in contact with the cooling medium, and “volume” is the
water at 81 8C; h = 420 J N s–1 N m–2 N K–1, kx = ky = 1.9 6 104
erg N s–1 N cm–1 N K–1 volume of the article.
According to its definition, k is 0.83 for a 3 mm 6 12
mm rectangular cross-section and 0.33 for a 6 mm thick-
ness slab. The higher value corresponds to a much more
cools down much slower and thus, resulting in the higher rapid cooling and moreover, it explains the much lower
crystallinity. This is the origin of the high gradient of crystallinity of the former compared to the latter.
crystallinity distribution. It is, however, not clear yet why Fig. 13 and 14 show crystallinity distributions of a 6
the crystallinity on the center is also very low. Moreover, mm 6 6 mm rectangular cross-section with k = 0.67 and
crystallinity can reach as high as 30% and 50% if a slab a 6 mm 6 24 mm rectangular cross-section with k = 0.42.
is cooled in water at 71 8C and 92 8C as demonstrated in Both samples were cooled in water. It is obvious that
Fig. 5. Crystallinity of the 3 mm 6 12 mm cross-section crystallinity is related to k: the larger the value of k, the
is quite different from that of the 6 mm thickness slab. lower the crystallinity obtained. Maximum crystallinity
According to the evolution of crystallinity distribution and average crystallinity at different k values are listed in
by cooling in silicon oil and water, we can conclude that Tab. 1.
a mild-cooling environment is favorable to obtain articles Influence of material properties: In previous sections,
with uniform crystallinity distribution and high crystalli- we focused on the influences of processing conditions
nity. and geometrical shape of the articles on the crystallinity
174 D. Yan, H. Jiang, H. Li

Tab. 1. Relationship between k and crystallinity for cooling in

water at 81 8C

Geometry k Maximum Average crystal-

in mm 6 mm in mm–1 crystallinity in % linity in %

3 6 12 0.83 4.3 2.8

666 0.67 11.0 7.0
6 6 24 0.42 40.0 24.6
66v 0.33 41.2 33.4

trary, the coefficient of thermal conductivity k has a sig-

nificant effect (ca. 15%).
Except for material properties, all the parameters
involved in the kinetic equation of crystallization are
empirical ones and can be obtained by fitting to experi-
mental data. Their accuracy depends on the experimental
Fig. 13. Crystallinity distribution of a 6 mm 6 6 mm sample
data. Variation of these parameters corresponds to the
cooled in water at 81 8C; for 60 s (needed for equilibrium crys- shift or deformation of the K(T) vs. T curve shown in
tallinity); h = 420 J N s–1 N m–2 N K–1, kx = ky = 1.9 6 104 Fig. 1. We further calculated the crystallinity distribution
erg N s–1 N cm–1 N K–1 of the sample cooled in silicon oil at 81 8C, and all condi-
tions are identical with those in Fig. 9 except that the
crystallization rate constant K (T) was increased by 10%.
The crystallinity distribution obtained (omitted here) is
almost identical with that in Fig. 9. On the whole, slight
variation in the parameters of crystallization kinetics has
only marginal effects on crystallinity distributions. In
order to get uniform crystallinity distributions, more
attention should be paid to processing conditions, such as
cooling temperature and cooling strength, and geometri-
cal shape of the articles.

Conclusion
In this work, we investigated the evolution of crystallinity
distributions on two-dimensional space. The crystallinity
obtained in nonisothermal crystallization processes
depends on the cooling rates. Crystallinity decreases with
Fig. 14. Crystallinity distribution of a 6 mm 6 24 mm sample increasing cooling rates, as demonstrated in Fig. 2. How-
cooled in water at 81 8C for 120 s. (needed for equilibrium crys- ever, the influence is negligible as long as the cooling
tallinity); h = 420 J N s–1 N m–2 N K–1, kx = ky = 1.9 6 104 rate is below a critical value (about 30 8C N min–1 for
erg N s–1 N cm–1 N K–1
PET), which makes it possible to realize high and uni-
form crystallinity under nonisothermal conditions. On
increasing the cooling rate above this critical value, a ser-
distribution. It has been clearly demonstrated that their ious gradient of crystallinity distribution can be observed.
effects are predominantly significant. Now we are going The crystallization behavior in mild and strong cooling
to discuss the influences of material properties and crys- media is quite different from each other. In the mild med-
tallization kinetics on crystallinity distributions. ium silicon oil, crystallinity is very high and the crystal-
Usually, material properties are functions of tempera- linity distribution is fairly uniform inside the rectangular
ture and crystallinity and in practice, we always use cross-section. During the initial period, crystallinity on
approximate values. In order to clarify the dependence of the surface is higher than in the internal part, and a bowl-
crystallinity on material properties, we increased every like pattern is formed. Both on the surface and in the cen-
material property to 110% and calculated the correspond- ter, crystallinity increases slowly with time and finally,
ing crystallinity. Calculations show that density q, speci- the crystallinity of the internal section exceeds that of the
fic heat capacity Cp, and latent heat of crystallization L surface and forms an inverted bowl. In the strong cooling
exert little influence on crystallinity (a 5%). On the con- medium water, crystallinity is very low and a serious gra-
FEM simulation of nonisothermal crystallization, 1 ... 175

4)
dient of crystallinity distribution exists. The crystallinity A. Y. Malkin, V. P. Beghishev, I. A. Keapin, Z. S. Andria-
on the surface reaches a very low equilibrium value in a nova, Polym. Eng. Sci. 24, 1396 (1984)
5) A. Y. Malkin, V. P. Beghishev, I. A. Keapin, Z. S. Andria-
short time and changes little afterwards. Changing the
nova, Polym. Eng. Sci. 24, 1402 (1984)
shape of the articles increases the crystallinity of the 6) K. Nakamura, T. Watanane, K. Katayama, T. Amano, J.
internal section, but crystallinity on the surface remains Appl. Polym. Sci. 16, 1077 (1972)
7)
almost unchanged, which causes the significant gradient K. Nakamura, K. Katayama, T. Amano, J. Appl. Polym. Sci.
of crystallinity distribution. 17, 1031 (1973)
8) T. Ozawa, Polymer 12, 150 (1971)
The surface to volume ratio of an article is an important 9) J. Berger, W. Schneider, Plast. Rubber Process. Appl. 6, 127
index, which can be used to estimate the temperature dis- (1986)
tribution and evaluate the crystallinity and its distribution. 10) T. W. Chan, A. I. Isayev, Polym. Eng. Sci. 34, 461 (1994)
11)
It can also be used to analyze the influence of the geome- G. Eder, H. Janeschitz-Kriegl, S. Liedauer, Prog. Polym. Sci.
trical shape as demonstrated in Tab. 1. 15, 629 (1990)
12) A. I. Isayev, T. W. Chan, K. Shimojo, M. Gmerek, J. Appl.
Polym. Sci. 55, 807 (1995)
13) A. I. Isayev, T. W. Chan, M. Gmerek, K. Shimojo, J. Appl.
Acknowledgement: This work was sponsored by the National Polym. Sci. 55, 821 (1995)
14) A. Mehl, L. Rebenfeld, J. Polym. Sci.: Part B: Polym. Phys.
Sciences Foundation of China and National Key Projects for
Fundamental Research “Macromolecular Condensed State” of 31, 1677 (1993)
15) A. Mehl, L. Rebenfeld, J. Polym. Sci.: Part B: Polym. Phys.
The State Sciences and Technology Administration of China.
31, 1687 (1993)
16) N. Billon, J. M. Haudin, Colloid Polym. Sci. 271, 343 (1993)
17) N. Billon, C. Magnet, J. M. Haudin, D. Lefevre, Colloid
1) M. Avrami, J. Chem. Phys. 7, 1103 (1939) Polym. Sci. 272, 633 (1994)
2) M. Avrami, J. Chem. Phys. 8, 212 (1940) 18) G. C. Alfonso, personal communication
3) M. Avrami, J. Chem. Phys. 9, 177 (1941) 19) H. Jiang, M. S. Thesis, Shanghai Jiao Tong University, 1998