. .... a... .
5...!

 

 

 

 

 

 

 

§t£fi1

Ely-:1.” t
x ’ ‘
l

.mu.

 

 

 

 

 

 

 

p .12.!
.Al) .
.L - I--.l.aul-n
90,2111: Linn-'1
i i. .5 l a. ‘
‘. z, riui .w

 

 

 

 

 

 

, .. ‘ 1O} 4
327.! .E ‘33:..qu .
.v Pz)x‘lm\\f . fix. .4. .1
V, i .tllztiirtt‘lli V”

qu..mmuEgZWdewru . y . , A . V. . .. . .. . ‘ , » ‘ . “ ‘ . , .1
9. ; fikfiefi, , , I ‘ , , ::. ._ . , ,. .1? 25:5".

 

ppm

WWIlllllllllllllllllllll!”ll"HI!HIIHIHIUHHIW

302074 1850

 

 

LIBRARY

Michigan State
University

 

 

PLACE IN RETURN Box to remove this checkout from your record.
To AVOID FINE return on or before date due.
MAY BE RECAU£D with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

11/00 WWW.“

MODELING THE SOLIDIFICATION OF
SEMICRYSTALLINE POLYMERS

By

Durgesh R696

A THESIS

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

MASTER OF SCIENCE

Mechanical Engineering

1998

ABSTRACT
MODELING THE SOLIDIFICATION OF
SEMICRYSTALLINE POLYMERS
By
Durgesh Rege

Modeling the solidification of semicrystalline polymers is considered in this work.
This is done to study the possibility of fabricating in-situ polymer/ polymer compos-
ites using directional solidification. With this technique, a thermal gradient is applied
to induce orientation of the crystals in a polymer blend composed of a semicrystalline
polymer and an amorphous polymer. A numerical model is first developed based
on the Eikonal equation to track the interfaces of the developing spherulites. The
travel times of the growing spherulites can be computed rapidly using a simple finite-
difference scheme. Isothermal solidification with instantaneous nucleation is consid-
ered for two-dimensional crystallization processes. The model correctly accounts for
the growth of single or multiple spherulites and their impingement without special
treatment. The model also accounts for the growth of spherulites in the presence
of obstacles such as fibers or the growth of spherulites in the presence of a thermal
gradient. Theoretical work is also performed to identify the thermal conditions that
favor steady growth and suppress excessive nucleation ahead of the growing inter-
face. Furthermore, stability issues are addressed to estimate whether a solid/ liquid

interface will be planar under directional solidification conditions.

To my family

iii

ACKNOWLEDGMENTS

I would like to express my deepest appreciation to the following people for having
made the completion of this thesis possible :

To my advisor Dr. Andre Benard, for his insight, guidance and invaluable assis-
tance, without which this thesis would not have been possible. His dedication to the
subject and personal integrity have served as an example and an inspiration to me.

To my other committee members, Dr. McGrath and Dr. Somerton for their
valuable suggestions and comments.

To Research Excellence Funds and the Composite Materials and Structures Center
for providing financial support for two years.

To all my other instructors at MSU, for their skillful instruction in the classroom,
but more importantly, for their constant encouragement and friendship.

To my friends, for having made my stay here the memorable experience that it
has been.

And to my family, for years of support, patience and understanding.

iv

TABLE OF CONTENTS

LIST OF FIGURES vii
LIST OF TABLES ix
1 Introduction 1

2 Crystallization in the Presence of a Temperature Gradient and So-

lute 9
2.1 Introduction .................................. 9
2.2 Nucleation of Spherulites ........................... 10
2.3 Growth of spherulites ............................. 15
2.4 Directional Solidification ........................... 18
2.4.1 Problem Formulation and Calculations .................. 21
2.4.2 Results and Discussion ........................... 23

3 Morphological Instability of a Planar Interface during Directional

Solidification 25
3.1 Introduction .................................. 25
3.2 Computation of the Solute Concentration Ahead of a Wavy Interface . . 27
3.3 Analysis of morphological instability .................... 31
4 Using the Eikonal equation to Model the Phase Change in Semicrys—

talline Polymers and their Composites 34
4.1 Introduction .................................. 34
4.2 Factors influencing the Spherulitic morphology ............... 35
4.3 Previous work done in numerical modeling ................. 38
4.4 Basis of the simulation ............................ 39
4.5 Vidale’s scheme ................................ 43
4.6 Sethian’s approach .............................. 47
4.7 Implementation of the scheme ........................ 50
4.8 Results and Discussion ............................ 51
4.8.1 Spherulitic Growth in a Neat Polymer .................. 51
4.8.2 Spherulitic Growth in the Presence of Obstacles ............. 55
4.8.3 Spherulitic Growth in the Presence of a Temperature Gradient ..... 58
4.9 Conclusions .................................. 58
5 Summary and Conclusions 61

v i

6 Suggestions for Future Work 63
APPENDICES 65
A C code based on Sethian’s approach 65
B Matlab program for simulations based on Vidale’s scheme 78
B.1 mex.m ..................................... 78
B.2 chmi.m ..................................... 81
B.3 middlem ................................... 87

B.4 corner.m .................................... 88

1.1

2.1
2.2
2.3

2.4
2.5

3.1
3.2

3.3

4.1

4.2

4.3
4.4

4.5
4.6
4.7
4.8

4.9

LIST OF FIGURES

Micrography of a spherulite taken in the laboratory of Dr. McGrath (Cour-
tesy of Eric Kaiser) ...........................

Nucleation rate as a function of temperature of crystallization ......
New nucleation rate as a function of temperature of crystallization

Growth rate (m/sec) as a function of temperature of crystallization(°C)
for isotactic polypropylene ........................

Schematic for the directional solidification apparatus used ........

Temperature gradient (°C/cm) as a function of velocity of growth
front(m/sec)for isotactic polyprOpylene .................

Solute concentration ahead of the planar solid/ liquid interface ......
Solute concentration ahead of the perturbed solid/ liquid interface (x nd y

distance is in m) .............................
Stability criterion for isotactic polypropylene ...............

Initial concept (now deemed erroneous) used in modeling of a spherulite
nucleating at N and growing around a circular obstacle ........

A schematic diagram of the shadow behind the circular obstacle with cen-
ter at M ..................................
Illustration of a propagating circle and the corresponding zero level set
A schematic depiction of a source point A (nucleus) and its neighbors in
the ring ..................................
A 2-D grid Showing the numerical calculation of travel times .......
An illustration of the sequence of solution of one edge of the ring
Flowchart for Sethian’s method .......................
Front position of a spherulite at different times while growing in an uniform
temperature field (Grid Size 2 100x100 microns and the times are in
microns/ sec) ...............................
Shape of two spherulites at different times illustrating impingement and
subsequent growth (Grid size = 100x 100 microns and the times are in
microns/sec) ...............................

4.10 Impingement of multiple spherulites with randomly located nuclei shown

at different instants (Grid size : 100x100 microns and the times are
in microns/sec) ..............................

vii

15
16

19
20

24

30

3O
33

36

37
41

54

viii

4.11 Growth of a spherulite in the presence of a rectangular obstacle Shown
at different times (Grid Size 2 100x100 microns and the times are in
microns/sec) ............................... 56
4.12 Spherulitic growth in presence of multiple obstacles shown at different
times (Grid Size 2 100x100 microns and the times are in microns/sec) 57
4.13 Spherulitic growth in the presence of a temperature gradient shown at

different times (Grid size = 100x 100 microns, times are in microns/ sec
and dT/dx = 104°K/m ) ........................ 59

LIST OF TABLES

2.1 Values of the parameters from data published in Prog. Polym. Sci . . . . 23

3.1 Assumed values of various parameters ................... 29

ix

Chapter 1

Introduction

Among the multitude of materials available today, there exists a class of materials
called polymers. Their range of behavior and molecular configurations are so wide
that further sub-divisions are inevitably made. In fact so different are their prop-
erties and processing requirements that they form the basis of an industry of their
own. Polymers, or plastics as they are commonly referred to when man-made, are
a relatively new class of materials. Although first developed toward the end of the
nineteenth century, synthetic polymers found major practical uses only after World
War II mainly due to cut-off of key natural raw materials. Since then the polymer
industry is advancing in technology and volume growth and today is one of the most

prosperous industries in the field of materials.

Polymers can be distinguished broadly as thermoplastics or thermosets. Thermo-
plastics are fully polymerised in their raw state and there is essentially no chemical
reaction involved in processing. Application of heat results in softening or melting
of these materials at which time the material flows and can be formed or molded
into the desired Shape. This cycle is reversible and allows the reuse of a polymer in
various products. Thermosets on the other hand are not fully polymerised in their

raw state. They contain a partially polymerised polymer and a substance called a

1

2
catalyst. Application of heat and pressure first causes softening and after several
minutes the catalyst breaks down and initiates a chemical reaction which completes

the polymerization or “curing” of the polymer.

In this work the class of polymers under consideration is thermoplastics. It will
actually be a study of thermoplastic polymer blends. The growth of this blend tech-
nology has significantly increased in the last two decades for the obvious reason that
customized property profiles are not always available from a pure polymer at the
desired cost. It is always more time and cost efficient to modify the existing polymer

product through blending and alloying.

The microstructure of thermoplastic polymer products is one of the dominant
factors that determine their chemical and mechanical properties. The properties of
semicrystalline polymers, in particular, are generally more sensitive than other materi-
als to the processing conditions as their microstructure can be influenced dramatically

by the processing method chosen.

The microstructure of these semicrystalline materials is generally defined by the
growth of crystalline entities called spherulites. The growth and morphology of
spherulites arising in semicrystalline polymers have been studied extensively over

the past thirty years. A brief literature review is outlined in the next section.

Some of the earliest work done on the microstructure of semicrystalline polymers
was by Keith and Padden [1, 2] at the Bell Telephone Laboratories. Their work gave
the first insight into the causes of the Spherulitic morphology observed in polymers.
In most examples of crystal growth, a primary nucleus grows into a crystallite having
a discrete crystallographic orientation. In spherulites however, primary nuclei also
initiate the formation of polycrystalline aggregates with no crystallographic orien-
tation, which are more or less radially symmetric. It was also discovered that all

high density polymers crystallize from the melt or from concentrated solution with a

 

Figure 1.1: Micrography of a spherulite taken in the laboratory of Dr. McGrath
(Courtesy of Eric Kaiser)

Spherulitic habit and hence there has been a significant interest in this field. It was
also noticed that there are similarities between Spherulitic crystallization in a wide
range of substances - organic and inorganic, polymeric and nonpolymeric etc. These
similarities were so pronounced that they could scarcely be coincidental. Hence a
conclusion was reached that all of the mechanisms of Spherulitic crystallization could

be accounted for on a unified basis.

Keith and Padden[1] concluded that Spherulitic crystallization is a characteristic
property of all systems in which two specific requirements are satisfied. First, con-
ditions should be favorable for the formation of crystals with a fibrous habit and
secondly that the fibers formed must be capable of noncrystallographic “small an-
gle branching”. Individual fibers branch and develop into sheaves followed by further
branching and ultimately resulting in Spherulitic growth. They made a detailed study
of all the significant physical factors which are responsible for this Spherulitic crystal-
lization. One of the first facts observed was that the growth of a spherulite at constant

radial rate represents a steady state process and not an equilibrium one. The next

4
most important observation was the vital role played by the impurities in promoting
a fibrous habit in Spherulitic crystallization. An impurity rich layer of thickness 6
= D/G, where D is the diffusion coefficient for impurity in the melt, and G is the
radial growth rate of a spherulite, which is formed in the liquid at the interface plays a
major role in inducing a fibrous habit in these melts. Regarding, “small angle branch-
ing”, screw dislocations or local temperature gradients were thought to be the cause
for fibers to branch but the noncrystallographic nature of this branching could not
be explained. This feature was attributed to sufficiently small diffusion coefficients
which ensure that values of 6 are small enough for branching to be predominantly
noncrystallographic. In fact it was inferred that the frequency of noncrystallographic
branching varies inversely with the magnitude of 6. This study was a milestone in the
progress of polymer engineering. To summarize the observations made by Keith and
Padden [1], it could be said that: a spherulite can be regarded as a particular type
of saturation dendrite which rejects impurities which get trapped between the fibers.
The melt is extremely viscous and hence the growth of fibrils is fast in comparison
with the rate at which impurities can diffuse away. Hence kinetics of radial growth

are linear and the process is controlled by nucleation rather than mass diffusion.

However this theory proposed by Keith and Padden[1, 2] was challenged by Bassett
et al.[3]. Their studies concluded that solutal effects are not the fundamental reason
for the Spherulitic morphology as proposed by Keith and Padden. They obserevd
certain spherultic microstructures without the ‘fibers’ or cellular habits. Also, the
solute layer of thickness 6 does not always correlate with lamellar width as discussed
earlier. Hence it was proved that kinetic effects and crystallization temperature also
played a major role in the formation of spherulites and solute diffusion was not the

only cause of the Spherulitic morphology.

Extensive analytical modeling has been carried out over the past fifty years to de-

.5
scribe the crystallization of polymers under isothermal and non-isothermal conditions.
The earliest of the models predicted the evolution of the crystallinity as a function of
temperature. This was followed by modeling the Spherulitic growth in polymers based
on the concepts of metallurgy. These models accounted for the heat and mass bal-
ance as well as the interfacial kinetics. Avrami[4] followed by many others developed
analytical models to address the important aspects such as nucleation. The theory of
growth of polymer crystals (the lamellae of a spherulite) was originally developed by
Hoffman and Lauritzen[5]. Further issues such as growth of multiple spherulites and
the inevitable impingement between them was studied by Avrami who introduced a
correction factor that related the volume of spherulites growing freely in an infinite
medium to the actual volume of impinging spherulites. However all these models are
limited in their applications. They cannot predict growth morphology in polymer
composites. They cannot account for phenomena such as growth of spherulites in the
presence of one or more obstacles or ‘fibers’ as well as in the presence of a thermal
gradient. Also, they fail to address issues at the microscopic level such as the stability

of the solid / liquid interface.

Landmark work was done by Mullins and Sekerka [6, 7] in 1963 when they estab-
lished the stability criterion for the solid-liquid interface. They theoretically investi-
gated the stability of the Shape of a moving planar liquid-solid interface during the
unidirectional freezing of a dilute binary alloy by calculating the time dependence of
the amplitude of a sinusoidal perturbation of infinitesimal amplitude introduced into
the planar shape. They deduced the stability criterion which was expressed in terms
of growth parameters. Instability occurs if any arbitrary perturbation grows and sta-
bility occurs if all components decay. Using parameters such as surface tension, the
conditions of phase of equilibria, temperature gradients on both Sides of the interface

and the solute concentration gradient, they derived the conditions of stability as a

6
function of the wave number of infinitesimal sinusoidal perturbations. Later Trivedi
and Kurz [8] extended the linear perturbation theory to the case of large thermal
Peclet numbers. Coriell et al. [9] studied the morphological stability of a binary alloy
during the initial transient period of directional solidification , in which the interface
velocity, concentration and temperature gradients are changing with time. By intro-
ducing sinusoidal perturbations of the planar crystal-melt interface they numerically
calculated the time evolution of these perturbations. Glicksman [10] reviewed the
basic concepts of heat and mass transport surrounding a growing dendrite and also
reviewed the dendritic scaling laws and influences of solutal and thermal diffusion
on dilute binary alloys. Recently, Schultz et al. [11] studied the incipient instability
of growth in polymer blends. They formulated a balance between the breakdown
of the solid-liquid interface and stabilization which determines the critical radius of
instability of a growing Sphere. Their instability analysis includes the Lauritzen and

Hoffman interface kinetics and capillarity effects as well.

In this work, the growth of spherulites is modeled first. Situations such as growth
in the presence of obstacles and in a temperature gradient are considered with a
numerical model. The analytical crystallization models discussed above cannot ac-
count for all of these aspects. Such studies have been performed in the past, but
the approach in this work is novel and very efficient. A brief outline of the previous
work done in modeling associated with polymer processing is presented. Lovinger
and Gryte [12] developed a numerical model for the shape of polymer spherulites
formed in a temperature gradient. Their model considered the polymer sample to be
divided into a series of narrow isothermal regions and, using an analogy to Snell’s law
[13], correctly predicted the observed fibrillar curvature. Their computer simulation
results were consistent with the experimental results. Similarly, Huang and Peter-

mann [14] performed some experiments to model the growth of oriented spherulites

7
of polybutene-l in a temperature gradient. When crystallized in a thermal gradient,
there is suppression of the nucleation process in the high temperature region and
the Spherulitic growth is obstructed except towards the higher temperature. They
employed a polarizing microsc0pe to study the development of spherulites and demon-
strated that a Single measurement of the spherulite growth in the temperature gra-
dient was sufficient to calculate many crystallization isotherms. More literature can
be found regarding the modeling of spherulites. In this work, modeling of spherulites
is carried out using a new concept called the level-set method. Level set methods are
numerical techniques which can follow the evolution of interfaces. Rather than follow
the interface itself, the level set approach takes a curve representing the interface and

builds it into a surface which is called the level set function.

In the second part of this work, a model is derived to account for the diffusion
of non-crystallizable Species in a blend subjected to a temperature gradient. This
is done to determine the possibility of creating in-Situ polymer/ polymer composites.
The Simple case of a linear temperature gradient is studied since it is possible to study
such systems in the laboratory. Previous research was done on the directional solid-
ification of polymers as it is a situation which is encountered in practical problems.
Diffusion of non-crystallizable components or segregation of impurities is almost in-
evitable. It was Keith and Padden [1, 2] who first systematically studied the influence
of impurity segregation on the kinetics of Spherulitic crystallization. Much later, Ryan
and Calvert [15] studied diffusion and annealing in crystallizing polymers followed by
which Billingham et al. [16, 17] made a detailed study of the diffusion of stabilizing
additives in polypropylene and also of atactic fractions in isotactic polypropylene.
The details of all these studies will be mentioned in the subsequent chapters. It is
hoped that the model which is developed in this work will allow the modeling of

microstructure evolution in more complex systems.

8
This work constitutes the analytical part of the work funded by Research Excel-
lence Funds. This thesis is divided into six chapters the first of which is this chapter
of introduction. The subsequent chapters will focus on the actual work done in this
project including some experimental observations related to the directional solidifica-

tion of the polymer samples.

Chapter 2

Crystallization in the Presence of a

Temperature Gradient and Solute

2.1 Introduction

Prediction of microstructure formation in polymers is of technological importance as
their properties depend on it. A number of studies have therefore been performed
to predict the final microstructure of a part and its temperature during cooling. An
approach based on the direct numerical Simulation of the growth process has been
used in the previous chapters. The boundary of each growing spherulite was followed
throughout the solidification process and numerical simulations were performed for
a small representative volume. In this chapter, the formulations for a theoretical
description of the Spherulitic growth process in the presence of a temperature gradient,
and in the presence of a solute are presented. This material is needed as it is one of
the objectives of this study to determine the conditions needed to suppress nucleation
ahead of the solid / liquid interface. It is hoped that a cellular morphology will result
from the instabilities created by the solute and by suppressing nucleation ahead of

the interface.

10

To begin the analysis, it is important to review the main stages of any non-
isothermal process of microstructure formation. These stages are:

1] Induction time period: In this time period, the polymer chains reorganize
themselves to form nuclei.

2] Nucleation: This is the dominant process at the beginning of solidification and
leads very rapidly to the establishment of the final spherulite density.

3] Growth: This is an important feature of polymer morphology Since it could
give rise to a variety of morphological units.

4] Impurity segregation: Impurities play a major role in determining the crystalline
morphology Since they are rejected preferentially and their diffusion determines the
overall morphology.

Each of these steps will be discussed below and will be supported by relevant
case studies. An analysis is then presented to compute the temperature gradient to

suppress nucleation ahead of the solid / liquid interface.

2.2 Nucleation of Spherulites

Spherulites are the largest microstructural features encountered during the solidifica-
tion of a quiescent semicrystalline polymer melt. They are polycrystalline aggregates
formed from a radiating array of crystalline “ribbons” that branch regularly to cre-
ate a three-dimensional structure of approximate radial symmetry. According to the
mechanism of Spherulitic crystallization, after nucleation, lamellae start growing ra-
dially to form axialites. These axialites which are the precursors of Spherulites, start
from stacks of lamellae with a planar geometry to eventually evolve into a full spher-
ical entity. Before delving much into spherulites, an indepth study of the nucleation
process will be presented which is responsible for their formation.

Nucleation is generally referred to as the formation of stable nuclei that are able

11
to grow for a given set of external conditions. Nucleation can be broadly classified
into four categories:

1] Homogeneous Nucleation: This is referred to as the Spontaneous formation of
nuclei from the aggregates of polymer chains, when the former exceed some critical
Size.

2] Heterogeneous Nucleation: Nucleation is said to be heterogeneous when nuclei
are formed on solid particles already existing in the polymer.

3] Thermal Nucleation: A process that occurs throughout the crystallization is
referred to as thermal nucleation.

4] Athermal Nucleation: This type of nucleation characterizes a solidification pro-
cess where all the crystals start growing at the same time.

Homogeneous nucleation is connected to thermal nucleation whereas heteroge-
neous nucleation can be either thermal or athermal. In short, the terms such as ho-
mogeneous and heterogeneous characterize the origin of the nuclei while terms such
as thermal and athermal characterize the temperature dependence of the nucleation
process.

The theoretical treatment of homogeneous nucleation occuring at irregular inter-
vals is based on the ideas developed by Turnbull and Fisher [18], and results in a
relation between nucleation rate, N, and temperature T of the following form:

—ED — AG”

N = No e$p[ RT

 

], (2.1)

where N, is a material constant and is proportional to temperature T, ED is the
activation energy for transport across the surface of the nucleus, AG" is the free-
energy of formation of the critical size nucleus and R is the universal gas constant.
For a normal three-dimensional nucleus AG“ is proportional to £17295, where Tm is the

melting-point, and AT is the degree of super-cooling. At small degrees of supercooling

12
(5 to 20°C) it is thus predicted that the nucleation rate will increase very rapidly
with decreasing temperature of crystallization, after an initial period of insensitivity.
At much lower temperatures, the EEE term dominates, so that N passes through a

maximum, and subsequently decreases with temperature.

However, homogeneous nucleation almost never occurs in polymers, because of
the fact that no polymer melt is pure. It always contains some form of impurity,
some nucleating agents or some such solid material. Hence the nuclei tend to form
on these existing solid particles thereby characterizing a heterogeneous nucleation.
Cormia et al. [19] have shown that heterogeneous nucleation is the dominant mode
of nucleation in most of the polymer systems. They used a very ingenious method
to support this evidence. Using a fine dispersion of tiny droplets of polyethylene,
they found that some particles crystallized at high temperatures of about 125°C , but
the majority which were free from heterogeneities or nucleating centers crystallized
at much lower temperatures of 90°C. Since the degree of undercooling needed for
homogeneous nucleation is far greater than for heterogeneous nucleation, the former

is not observed in the practical crystallization of the bulk material

Although heterogeneous nucleation can be either thermal or athermal, the latter
is most predominant and will be influencing the nucleation process under study. The
modeling of heterogeneous nucleation is a complicated process, since it is necessary to
know the distribution of potential nuclei in advance which depends on a large number
of factors such as the crystallization temperature, the processing conditions etc. The
process is further complicated by the fact that a distinction is made between the
apparent and true nucleation rate. While apparent nucleation is found from direct
observation of the Specimen, the latter is the rate of nucleation of untransformed
material accounting for those nuclei which are engulfed by the growing Spherulites as

well. However, to avoid further complications, it is assumed that the true nucleation

13

rate is equal to the apparent nucleation rate.

The fact that spherulites in most polymer systems undergo heterogeneous nucle-
ation was further substantiated by Chew et al [20]. They found that the Spherulites
reappeared in almost identical positions when the sample had been melted and re-
crystallized. Secondly, both the nucleation rate measurements and the Spherulite Size
distribution diagram confirm that the majority of nuclei form in early stages of trans-
formation which implies that the degree of undercooling is considerably less. During
these early stages of transformation, the number of nuclei increases almost linearly
with time. It has already been Shown by Cormia et al [19] that the nucleation occur-
ing at low degrees of undercooling is essentially due to the presence of heterogeneities

and is described as heterogeneous or pseudo-homogeneous nucleation.

Nucleation also can be categorized as isothermal and nonisothermal, depending
on the processing conditions. In isothermal nucleation, according to Avrami [4], there
is a uniform and fixed density, M, of potential nuclei and they nucleate randomly in
time with each nucleus having a probability, V, of developing per unit time. This can

be written as

N(7') = u.Me:rp(—VT), (2.2)

where N is the rate of nucleation and T is the time function. The nucleation rate is
thus an exponentially decreasing function of time. The total number of nuclei as a

function of time can thus be written as

N(T) = M[1 — 8$p(—V7’)]. (2.3)

However, Icenogle [21] found that this equation did not match the experimental data
and hence suggested a relation depending on the time of onset of nucleation and

time of cessation of nucleation. Many more theories have been developed to describe

14
isothermal nucleation, though all of them require a Significant amount of experimental

data.

The nucleation process under study is non-isothermal and is most commonly en-
countered in polymer systems. A model is presented here to describe the nonisother-
mal nucleation process in which a temperature gradientis imposed on the system. Nu-
cleation can be predicted to a reasonable accuracy with semi-empirical rules, provided
all the data for a specific polymer is known, viz. nucleation rate, density, parameters
of the nucleation law used, etc. The phenomenon of nucleation depends mainly on the
thermal fluctuations which lead to the creation of variously sized spherulites and it
also depends on the creation of an interface between the liquid and the solid. Thus it
is the thermal fluctuations which lead to the formation of Spherulites. The nucleation

mechanism is represented by a temperature dependent function given by

N = N, exp[—[3N(Tc — T,)], (2.4)

where N, and [3N are fitting parameters. N, and. T, represent the reference nucleation
rate and the reference temperature respectively and T C is the actual temperature of
crystallization. This expression was found to give a satisfactory fit to experimen-
tal data describing Spherulitic nucleation in polypropylene blends[22]. The result

obtained by using published data[22] is shown in Figure 2.1.

The values for N, and flN are 2 x 10’1 /m3 and 0.18/°K respectively. The reference
temperature value was set to 115°C. AS per the figure, the nucleation rate is seen
to be tremendously high for the corresponding range of crystallization temperatures
for the polypropylene samples. Such a high rate of nucleation can be attributed to
the fact that any commercial polypropylene sample contains one or more nucleating
agents which results in the exceptionally high rate of nuclei formation per unit volume

of the sample.

 

1e+10

I

8e+09

I

6e+09

I

4e+09

N (nuclei/mAS)

I

Ze+09

 

 

 

l l 1 l

o x
130 135 140 145 150 155 160

 

To in C
Figure 2.1: Nucleation rate as a function of temperature of crystallization
Since the polypropylene sample under study is devoid of any kind of nucleating
agents, the reference rate of nucleation was reduced by two and the nucleation rate
was recalculated. The results so obtained are Shown in Figure 2.2 which seem to be

reasonable

2.3 Growth of spherulites

One of the most important features of polymer morphology is the growth process

which can give rise to a variety of morphological units. In devising the models of

growth mechanisms for polymer crystallization, it is necessary to understand dif-

ferent growth mechanisms which happen to be quite complex in nature in order to

establish the most dominant among them for a given set of conditions. During crys-
tallization from the melt, growth may not only involve branching fibrils with defects
and amorphous regions, but also it may pass through more than one stage during the
course of the process.

The mechanism of Spherulitic growth was described as early as 1929 by

16
1e+08

 

89+07

I

69+07

I

I

4e+07

N (nuclei/mAS)

29+07 r

 

 

 

 

o l l I l
130 135 140 145 150 155 160

Tc in C
Figure 2.2: New nucleation rate as a function of temperature of crystallization

Bernauer[23]. According to this mechanism, after nucleation, lamellae start grow-
ing radially to form axialites, which are precursors of spherulites. The axialites grow
to eventually form Spherulites by a process in which they start from stacks of lamellae,
with a planar geometry, to evolve into a full spherical entity. The work done by Keith
and Padden[1] provided the first phenomenological theory of Spherulitic growth and
experimental observations. They studied the growth of spherulites from polypropy-
lene melts. As discussed before, they postulated that two conditions were required
for the formation of spherulites: high viscosity of the medium and the presence of
non—crystallizable material. In an isotactic polymer with atactic chains, the atactic
portions crystallize leSS readily than the isotactic components and are rejected from
the growing solid. This rejected component is commonly refered to as the “solute”.
This solute rejected from the crystal creates an excess concentration of impurity which
is pushed ahead of the solid / liquid interface. This layer of solute has a bearing on the
Spherulitic growth. It also has the dual effect of locally depressing the equilibrium

crystallization point and reducing the growth rate due to a shortage of crystallizable

17
material. Thus modeling the growth of spherulites requires that the effects of nu-
cleation, interfacial kinetics, temperature and solute are all accounted for. A brief

description of the growth process is given below in the next section.

After the appearance of the nuclei, the growth process is initiated. The growth
velocity of Spherulites is dominated by interfacial kinetics, rather than by a heat
difl'usion process [24]. The radial growth velocity in addition to the crystallization
temperature, Tc, is also a function of the composition of the melt in the case of
a multi-component system. In case of pure, low molecular weight substances, the
growth velocity of spherulites may have a significant magnitude while N, the rate of
nucleation, may be very small. This region is referred to as the metastable zone of
low subcooling, in which the growth of existing nuclei may proceed without further
homogeneous nucleation. The mechanism of growth for Spherulitic bodies involves
the formation of secondary nuclei. The most important features of any Spherulitic
growth are the presence of a high viscosity medium and the presence of some dis-
solved impurity which can be rejected from the growing front. These impurities
could be added deliberately or unintentionally present. In polymers though, they are
present as branched, atactic or low-molecular weight components, which crystallize
less rapidly than the bulk polymer and hence are capable of being rejected during
growth. The original theory describing the kinetics of growth was developed by Hoff—
man and Lauritzen[24]. The growth rate of each branch can be represented in a
general manner by the equation

8R -—U* —Kg

7%— : Go exp[————R(T _ Too)] exp[TAT ]. (2.5)

 

The first exponential term represents the temperature dependence of the segmental
jump rate in the polymer, while the second term is the contribution of the net rate of

attachment of new chains to the crystal. R is the spherulite radius and Go is a pre-

18
exponential constant that combines terms not strongly dependent on the material.
The growth rate of a polymer crystal therefore results from the rate of deposition
of nuclei of a stable size and their subsequent growth, as well as the rate at which
untransformed chains are brought to the growing crystal faces.

Similar to the nucleation model, an analytical model is developed for the growth
of Spherulites. It has been proved by earlier work done that the growth process is
controlled by a secondary nucleation step. This can be deduced from the temperature-
dependence of the observed growth rates at small degrees of super-cooling. Regardless
of the nature of the nucleation process, the growth rate of a polymer spherulite can

be described by the equation

G 2 G, exp[—flG(Tc —— T,)], (2.6)

where G, and fig are fitting parameters. C, and T, represent the reference growth
rate and the reference temperature respectively and TC is the temperature of crystal—
lization. This expression was found to give a satisfactory fit to experimental data on
the rate of Spherulitic growth in polypropylene blends[22]. The result obtained by

using some published data[22] is Shown in Figure 2.3.

2.4 Directional Solidification

Many attempts have been made to obtain a Spherulitic morphology with specific ori-
entation. These techniques involved application of a Shear stress on the crystallizing
melt and a restriction of the nucleation rate of new spherulites in the melt. Price et
al.[25] crystallized PEO from a single nucleation point within a very thin capillary.
Sasaguri et al.[26] unidirectionally crystallized a thin film of polybutene—l in a tem-

perature gradient. Fujiwara et al.[27] and Crissman [28] independently crytallized

19

 

 

 

 

1e-05 . I 1 1 a
9e-06 ~ a
8e-06 - -
79-06 - -
0) 69-06 - -
259-06 - -
(5 49-06 - a
3e-06 - ]
26-06 - -
1e-06 - -

0 J 1 1 1 n

110 115 120 125 130 135 140
Tc in C

Figure 2.3: Growth rate (m/sec) as a function of temperature of crystallization(°C)
for isotactic polypropylene

polypropylene and polyethylene samples using the technique of zone solidification.
The process employed in this work is called directional solidification. It has been
used for several years in metals and this technique can be adapted to polymer blends
in order to manufacture in-Situ polymer/ polymer composites with enhanced proper-
ties compared to the original constituents. It was used by Lovinger et al. [12] to study
crystallization of polymers, but no attempt was made to manipulate the microstruc-
ture. Directional solidification consists in using a stable linear temperature gradient
in such a way, that one part of the sample is above the phase change temperature and
the other below it. When the sample is moved through the temperature gradient at
a constant velocity the solid / liquid interface remains at the same position and can be

observed continuously. This principle is illustrated by a schematic as shown in Figure

2.4.

20

INTERFACE
COVER SLIDE

/

THOT (T h )

 

I

I

>1
7
I

I

 

I'l--'|'l|'|lu|l" ’I'II'IIII'I'I'l-II'I

 

I
14
I
I
I
I
I
I
I
I
I
I
I
I
I
I
I
I
I
I
I
l
I
I
I
I
I
I
I
I
I
I

 

TCOLD (T c)

 

SAMPLE

 

 

 

 

(sample is pushed to the left)

'IIllrll'

 

 

 

 

 

 

Figure 2.4: Schematic for the directional solidification apparatus used

21

2.4.1 Problem Formulation and Calculations

The key to obtaining an uninterrupted growth front from a single nucleation point
lies in the suppression of the nucleation within the subcooled melt ahead of the grow-
ing crystal front. To suppress heterogeneous nucleation, low degrees of undercooling
must be used. However, there is a limit to the degree of undercooling that can be
employed while at the same time having a reasonable rate of crystallization. In this
work, an analysis has been done to determine the case for which Single crystal growth
can be maintained under conditions where the nucleation rate is not completely sup-
pressed. The analysis is a modified version of the work done by Lovinger et al.[12].
When a molten polymeric material is traversed through a linear temperature profile
that brackets its crystallization temperature, crystal growth will be oriented in the
direction of the temperature gradient, resulting in long, oblong Spherulites. The ther-
mal gradient minimizes nucleation in the high temperature region and this results
in oriented spherulites during the crystallization process. A polymer melt sample at
temperature, Tm, and having a cross sectional area, A, is considered for the analysis.
If the sample is moving at some velocity, u, through the heated zone, then at a tem-
perature, Tc, the growth rate, G, is identical with the velocity, u, of the sample. If it
is assumed that the temperature gradient, Ta, is constant over the region of interest,

then the temperature at any position, x, below the melting zone is given as:

Tc 2 Tm + T0515. (2.7)

This is the temperature of crystallization where the melt will start crystallizing to

form Spherulites.

The total number of primary nuclei formed, 1%, in a sample of cross-sectional

22

area, A, as it moves from melting point to crystallization point at a velocity, u, is

n,==NI: (28)

where V is the volume of crystallized polymer and N is the total number of nuclei

given by equation (2.4). This equation can also be written as

IIT = A/ Ndx, (2.9)
0
which is the same as,
A
m._ ——/T NdT. aim
Ta ,,,

Equation (2.10) is obtained due to a change of variables in equation (2.7). From
equation (2.7), it is clear that :I: = 19%;“? The variables are changed from x to Tc. N
could be substituted by the expression obtained in equation (2.4). Substituting the

value of N and the limits of integration, the equation obtained is

_.N,41

nr 2 TG (3N [exp —’3N (T T, )) — exp(——,I3N(Tm —- T,))] (2.11)

 

The critical condition for a Single crystal to be formed that is continuous and devoid

of subsequent spherulites is that nT be less than unity. Therefore,

nT < 1 (2.12)

From equation (2.6),
(2.13)

23

Table 2.1: Values of the parameters from data published in Prog. Polym. Sci

 

 

 

 

 

N, ‘2 x 105n—3
1‘30 0.916K—1

3N 0.18K-I

C, 1.259 x 10-6 m/s
T, 115°C

 

 

 

 

From equations (2.11) and (2.13), equation (2.12) can be written as

 

—AN, G {343.
131v [(G—) "G — explfi~(Tr — Tm)] < Ta. (2.14)

To compute equation (2.14) values of all the constants and variables that appear
in the equation are needed. The values Shown in Table 2.1 are known from data
published in [22].

In addition, the cross-sectional area of the sample is assumed to be 1 x 10'°m2,
the melting point as 200°C and the growth rate which is same as the velocity of the
sample as 10‘6 m / sec. Substituting all of these values into equation (2.14) the temper-
ature gradient Should be approximately 110°C/cm, as seen in Figure 3.5. From this
equation it is clear that at small cross-sectional areas the volume of the crystallized
polymer will be considerably reduced. Hence the solidification carried out in small

capillaries or between two narrow glass plates will lead to uniaxial crystallization.

2.4.2 Results and Discussion

It was observed experimentally that the spherical symmetry of the spherulites is
distorted in the direction of growth. AS the growth rate was reduced, the elongation

became more pronounced and the number of nucleation points progressively reduced.

 

 

 

 

'60 T r I l I r l T I
-30 .. -
E -100 ~ ~
9
O
.5
o -120 r' -*
l_|
-140 — -
-160 — -
l l J l 1 l L l l
0 19-06 29-06 39-06 49-06 59-06 69-06 79-06 89-06 99-06 19-05
G in III/sec

Figure 2.5: Temperature gradient (°C/cm) as a function of velocity of growth
front(m/sec)for isotactic polypropylene

At very low velocities, spherulites of unlimited length and oriented in the direction
of solidification were grown from the initial nuclei. There exists Sporadic nucleation
of spherulites within the bulk of the oriented material. This is caused by persistent
heterogeneous nuclei. However, only very few of these serve as heterogeneous nuclei
for additional nucleation and nucleation is highly suppressed because of the thermal
gradient applied. An uniform thermal gradient is assumed across the thickness of the
sample. This was numerically verified by calculating the Biot number (hL/ k), where
h is the convection heat transfer coefficient, L is the thickness of the sample and k is
the thermal conductivity of the polymer sample. The Biot number (Bi) provides a
measure of the temperature drop in the solid relative to the temperature difference
between between the surface and the fluid. For Bi <0.1, it is reasonable to assume a
uniform temperature distribution across a soild at any time. It is assumed that h:5
W/m2K, L = 0.1 mm and k = 0.027 W/m K. For these assumed values, Bi 2 0.018.
Hence the assumption of a uniform temperature distribution is reasonable, but at the

limit of being acceptable. This issue can be addressed in future work.

Chapter 3

Morphological Instability of a
Planar Interface during Directional

Solidification

3.1 Introduction

During solidification of an alloy, there is interaction of heat and mass transfer pro-
cesses which can lead to an uneven solidification front. This phenomenon is responsi-
ble for the characteristic dendritic structure of alloys. The fundamental explanation
for the onset of uneven growth has been provided by the classical Mullins-Sekerka
theory of morphological instability[6][7]. This analysis is the key to our understand-
ing of morphological features of solid-liquid interfaces. They developed a rigorous
theory of stability by calculating the time dependence of the amplitude of a sinu-
soidal perturbation of infinitesimal initial amplitude introduced into the shape of the
plane. These perturbations are analyzed in terms of infinitesimal Sinusoidal Fourier
components. They postulated the fact that the interface becomes unstable if any

Sinusoidal wave grows and is stable if none grow. This theory accurately predicts the

25

26
absolute stability of a planar interface at high velocities under constrained growth
conditions. However, it suffers from a principal drawback. Their results predicted
that a planar interface would be unstable at all velocities in an undercooled melt.
This is not true for the case of rapid solidification. Experimental observations have
Shown the existence of absolute stability of a planar interface in undercooled melts
which have been rapidly quenched. Hence some attempts were made to extend their
theory to general cases where rapid solidification is encountered. Trivedi et al.[8] re—
examined the conventional constitutional supercooling criterion and proposed a more
general planar interface stability criterion which is valid for low as well as high growth
rate conditions. They derived conditions for absolute stability of pure and alloy un-
dercooled melts and applied their results to dendritic growth from pure undercooled

melt.

Among the recent work done on this topic, is the numerical model developed by
Coriell et al.[9]. They modeled directional solidification of a binary alloy during the
initial transient period in which the interface velocity, concentration, and temperature
gradients are changing with time. They introduced sinusoidal perturbations of the
planar crystal-melt interface and numerically calculated the time evolution of these
perturbations. The results obtained by them were in good agreement with the Mullins
and Sekerka analysis, provided instantaneous values of the temperature gradient and
solidification velocity are used in the analysis. Schultz et al. [11] examined the
results of morphological stability of a polymer crystallizing from a binary melt for
large undercoolings. They showed that, for a polymer crystallizing, with its local
interface kinetics and rejection of solute, there exists a transition from a compact
Sphere to a more Open one at the point of morphological instability. Jung et al.[29]
non-dimensionalised the results obtained by Mullins and Sekerka’s stability analysis.

By means of a Simplified asymptotic analysis of the non-dimensionalised growth rate

27
equation they derived a single expression for the wavelength of maximum instability.

Their work provided a concise presentation of the stability theory.

3.2 Computation of the Solute Concentration

Ahead of a Wavy Interface

Under steady-state conditions, the concentrations and solid / liquid interface morphol-
ogy are independent of time. This implies that the solution to the basic diffusion
problem is indeterminate and a whole range of morphologies is permissible from a
mathematical point of view. In order to distinguish the solution which is the most
likely to correspond to reality, it is necessary to find some additional criterion. Ex-
amination of the stability of a Slightly perturbed growth form is probably the most
reasonable manner in which to treat this Situation. The general features of a solid-
ification problem can be described as follows: a solid/liquid interface whose form is
defined by a given mathematical function containing one or more variable parameters
assumed to be advancing without change into the melt. AS it advances, heat and/or
solute are evolved at each point of the interface and diffuse into the solid and the
melt. The diffusing solute will build up ahead of the interface when k < 1, where k is
the solute partition coefficient (ratio of interface concentration in the liquid to that
in the solid) and form a boundary layer, while a uniform level of the solute, Co, is
supposed to exist at a sufficiently large distance from the interface. The boundary
layer can be characterised by the ratio of the diffusion coefficient to the growth rate.
The diffusion coefficient can be either the thermal diffusion coefficient (in case of the
rejection of heat) or the mutual mass diffusion coefficient (in case of the rejection of
solute). The thickness of the thermal boundary layer due to the rejection of latent

heat is large compared to the scale of microstructures that develop in the polymer

28
blend, which are normally in the micron and sub—micron range. Thus it can be safely
assumed that thermal fields have no role to play in the pattern formation, at least
for the case of small velocities of growth. The other boundary layer is due to the
solute diffusion, in this case the diffusion lengths are small and can scale with the
microstructures observed in polymer blends. Hence attention will be restricted here

only to the solute diffusion occurring in the liquid.

The basic equation governing any diffusion process is considered [5] and after
an extensive mathematical treatment, namely, applying the method of separation of
variables and then a suitable set of boundary conditions, a complete solution for the
solute distribution is obtained. This equation which represents the solute distribution
ahead of a planar solid/ liquid interface, advancing under steady-state conditions is
written as

Q _ V2

C = C0 +( k C0)ea:p[——D—], (3.1)

If the geometry of the problem is such that the interface form corresponds closely
to some simple shape, any exact solution which is available for the Simple Shape
can be assumed to be Similar to that for the Slightly different morphology. The
planar solid / liquid interface is then assumed to be changed and the solute distribution
calculations for this slightly perturbed solid / liquid interface are made which result in

an equation of the form

kVCCe Sin[wy]
Vp — Db

 

C = C0 + (92 — C0)e:rp[—%] +( )errp[—bz]. (3.2)

k

Comparing the above two equations, it is seen that the latter has the addition of a
perturbation term. In these equations, Co represents the initial concentration of the
solute, k is the solute partition coefficient, D is the diffusion coefficient in the liquid

and V is the velocity of the growth front. In addition, the parameters 0,, z, b and p

29

Table 3.1: Assumed values of various parameters

 

 

 

 

 

 

D 10‘12m2/sec
k 0.3

V 10‘°m/sec
Co 0.2

A 10’6 m

c 10‘6 m

 

 

 

 

are defined by the following equations:
GC 2 —V/D(Co/k — C0),

2 = csin(wy), where e is the amplitude of perturbation and w = 27I/A is the wave

number,
__ v v 2 .
b— §T+VE +w2 and
p=1-k.

Certain values have been assumed for some of the parameters in the above equa-
tions. It can be seen that after substituting all these values in the above two equations
the solute concentration C is obtained as a function of distance. The plots for the
variation of the solute concentration as a function of y for both the planar as well as

the Slightly perturbed interfaces are shown in Figure 3.1 and 3.2 respectively.

A parameter of major significance in this analysis is the quantity 6 = D/ G, where
D is the diffusion coefficient for impurity in the melt, and G is the radial growth rate
of a Spherulite. The quantity (5, whose dimension is that of length determines the
lateral dimension of the fibers. From Figure 3.2 it is clear that the lateral dimension
of 19-06 In between the tips of the dendrite correlates well with the actual dendritic

spacing seen in the Spherulitic growth.

30

 

C at unperturbed interface

 

 

 

035 A 1 l l l J 1 1 l

 

0 19-07 29-07 39-07 49-07 59-07 69-07 79-07 89-07 99-07 19-06
Distance in m

Figure 3.1: Solute concentration ahead of the planar solid/ liquid interface

 

 

 

 

 

I I I I = 5 5 1906
I s I a' i I i
I E I i = i =
'I ' i E i i i
’I '=. 'I i 5. i i ‘ 56"”
t . ‘. I . I i
I I 5 i I i :
l E .9 i E I I
I. .' l i I i I
y I I! ' : i I _
'1' i i i I i i 0 Y
i s s E i i
it :1 i i- i i i
'1, i i E 5 E i - -59-07
= i i I i 3
i i i . 5 E I E
I’ 5 I i i i i
I I 2 I I 1 I l l I 1 L1 I - l -1e-06

0 19-07 29-0739-07 49-07 59-07 69-0779-07 89-07 99-07 19-06

X

Figure 3.2: Solute concentration ahead of the perturbed solid/liquid interface (x nd
y distance is in m)

31
3.3 Analysis of morphological instability

Classical thermodynamic definitions of stability are inapplicable to the determination
of the morphology of a growing interface. Hence it is necessary to use heurestically-
based stability criteria. The Simplest assumption made is that the morphology which
appears is the one which has the maximum growth rate or minimum undercooling.
Depending upon the temperature gradient in the liquid melt at the solid/liquid in-
terface there may, or may not, exist a zone of constitutional undercooling. This zone
is defined as that volume of the melt ahead of the interface within which the actual
temperature is lower than than the local equilibrium solidification temperature. For
the existence of such a constitutionally undercooled zone, the temperature gradient
at the interface in the liquid Should be lower than the gradient of the liquidus tem-
perature in the melt which is simply the product of concentration gradient and the
liquidus slope, m. The expression for the concentration gradient is given by equation

(3.1). The equation governing this criterion can thus be written as:

_ mVACo

T .
G< D I (33)

where Ta is the temperature gradient at the interface in the liquid, m is the liquidus
Slope, V is the velocity of the growth front, AC0 = 99%;” is the concentration
difference between the liquid and solid solute contents, where Co is the composition
of the polymer blend and k is the solute partition or the distribution coefficient. D

is the diffusion coefficient in the liquid. This equation is satisfied with the values of

parameters assumed in Table 3.1.

This constitutional undercooling criterion provides a useful means for estimat—
ing whether a solid/liquid interface will be planar under directional solidification

conditions. This is also called the pseudo-thermodynamic criterion for the stability

32
analysis. The interface will always become unstable so long as this equation is being
satisfied. However, this method has several faults. Firstly, it ignores the effect of the
surface tension of the interface which will tend to inhibit the formation of perturba-
tions. Secondly, it takes into account only the temperature gradient in the liquid.
Thirdly, it does not give any indication of the scale of the perturbations which will
develop if an interface becomes unstable. Hence, it is necessary to suppose that the
interface has already been slightly disturbed and to see whether the disturbance will
grow or disappear. This method is the most famous Mullins-Sekerka criterion for
evaluating the interface stability. The mathematical treatment for this theory is very
intensive and has been extensively covered in the book by Kurz and Fisher[30]. The
mathematical treatment has not been presented here but the final result has been
presented for further analysis. The governing equation which is used for the analysis
is given as
5' V V

c : (mGC)(b — 5)(_W2F — TC + mGC), (3’4)

 

where F is the Gibbs-Thomson coefficient, ch is the temperature gradient in the
liquid and 6 2 fig. The parameterf, describes the relative rate of development of the
amplitude of a small sinusoidal perturbation. For very small values of w, the value
of this parameter is negative due to the curvature damping and the perturbation
will tend to disappear. For wave numbers greater than A,- and above, the sinusoidal
Shape will become more accentuated. The wavelength having the highest rate of
development is likely to become dominant. The reason for the tendency to stability
at high A-values is the difficulty of diffusional mass transfer over large distances. Thus
the Sign of the parameter determines the stability of the interface. If it is positive for
any value of w, then the perturbations with that wave-number will be amplified and
they tend to disappear if the value of the parameter iS negative. From the calculations

0

it was found that f > O. This has been graphically shown in Figur93.3. Hence the

33

 

 

 

 

 

C
.9 _
E
e
3 _
t .
0)
3 0.4 - ~
0
g 0.2 - -
S
2 0 _ ............................................................................................................................ _
(D
>
8 oz — —
B 4
a) - _ A
*0 ° F
g -0.6 - -
«.6
m -0.8 - -
CE

_1 l I 1 l

0 200000 400000 600000 800000 1e+06

Wave number w

Figure 3.3: Stability criterion for isotactic polypropylene

interface is unstable and the wave number is amplified. The analytical results obtained
indicate that the interface concentration at the onset of morphological instability are
accurately predicted by the classic Mullins and Sekerka analysis. Once the values of
the interface velocity and the temperature gradient at the interface are known, the
critical interface concentration can be calculated from the Mullins Sekerka analysis

and it is easy to find out if the interface would stabilise or become further unstable.

Chapter 4

Using the Eikonal equation to
Model the Phase Change in
Semicrystalline Polymers and their

Composites

4. 1 Introduction

Modeling the crystallization of polymers has significant industrial importance Since
the properties of a product depend on the microstructure developed during processing.
For a neat polymer, the main microstructural features observed under the microscope
are Spherulites. They are aggregates of crystalline plate-like lamellae which grow
radially outward from a nucleus.

The Spherulitic microstructure defines numerous physical and chemical properties
of polymers and its development during processing has been the object of intense
research for the past forty years. One of the first quantitative models developed

to describe the isothermal evolution of crystallinity in polymers was independently

34

35
derived by Avrami [4], Johnson and Mehl [31] and Kolmogorov [32]. The classical
analysis of Avrami is often successful in modeling the crystallization process but it
also has a number of limitations which have been reported extensively in the past [33].
Among these limitations are variations in the linear growth rate of the spherulites,
which may not remain constant at all times. Furthermore, a two-stage crystallization
can occur which invalidates the model, and finally the presence of reinforcement such
as fibers physically constrains the growth of Spherulites and modifies their shape.
Analytical models have been developed to predict the effect of fibers on crystallization
but they do not account for the “Shadow” (merging of fronts behind the obstacle)
created by the inclusions in the melt [34] [35]. There was therefore a need to develop
a model which could accurately predict the influence of fibers on crystallization and

explain many of the unusual features of these systems.

Computer simulation of polymer crystallization is an effective means for studying
the crystallization kinetics and morphology development of semicrystalline polymers
as they allow flexibility in nucleation and growth models. A Significant amount of
work has been carried out over the past twenty years and many numerical models
have been developed [36] [37] [38][39] [40]. The approach presented below is based on a

rapid technique for solving the Eikonal equation[41].

4.2 Factors influencing the Spherulitic morphology

The process of Spherulitic growth in polymers can give rise to a variety of morpholog-
ical entities. Extensive experimental and modeling work has been done in the past
to study the growth of crystalline forms in many polymer systems. From this work
there are three different aspects that affect the morphology: growth in an isothermal

field, growth in the presence of fibers and growth of Spherulites in the presence of a

36

 

Figure 4.1: Initial concept (now deemed erroneous) used in modeling of a Spherulite
nucleating at N and growing around a circular obstacle
temperature gradient. These aspects will be discussed briefly below.

Spherulitic growth in isothermal conditions has long been studied. For example,
the work done by Keith and Padden [2] provided the first phenomenological theory
of Spherulitic growth and experimental observations. The Spherulitic growth process
is initiated after the appearance of the nuclei and the growth velocity of spherulites
is often dominated by interfacial kinetics. The growth rate in any isothermal field
therefore results from the rate of deposition of secondary nuclei of a stable size and
their subsequent growth. A roughly circular entity is formed that may impinge with
growing spherulites until solidification is completed.

The growth of spherulites in reinforced polymers has also been an area of study for
quite some time. In the isothermal case, each spherulite possesses the same constant
growth rate. However, once a spherulite reaches an inclusion in the melt, its growth
behavior changes. Initially, the growth of the spherulite as depicted in Figure 4.1

was considered to be the correct description. It was assumed that the spherulite

37

T Upper shadow

  

 

<
Nucleus
N

Lower shadow

Figure 4.2: A schematic diagram of the shadow behind the circular obstacle with
center at M

would grow circularly, whereby the growth through the obstacle was allowed. The
growth fronts were always parts of a circle and after the obstacle was totally engulfed,
the shape of the growth front was not disturbed just as if no obstacle was present.
But the experimental evidence of Schulze et al. [42] as discussed in the previous
section indicates that the growth as per Figure 4.1 is incorrect. At all times, the
surface of the spherulite stays normal to the surface of the obstacle and it does not
remain spherical in shape once an obstacle is encountered. There are certain distinct
differences between the growth morphology of spherulites in a neat polymer and those
in a polymer composite which have been discussed in detail by Schulze and Biermann
[42]. These authors conducted controlled experiments for the case of one spherulite
growing around a circular obstacle which is shown schematically in Figure 4.2. It
was found that the growth fronts remain perpendicular to the obstacle beginning at
the points T in Figure 4.2. The growth is correctly described in the regions outside
the “Shadow” of Figure 4.2. The “Shadow” is limited by two straight lines which are

tangential to the border of the obstacle and which pass through the nucleus, N. The

38
tangent points are denoted by points T. The Shape of the spherulite is circular only
until it impinges with the obstacle. Once the tangential lamellae come in contact with
the obstacle, they grow into the liquid zone behind the obstacle and wrap around it.
The growth of spherulites in a thermal gradient was studied by Lovinger et al.[43],
among others. They conducted many experiments on samples of polymer systems by
directionally-solidifying them in temperature gradients up to 300°C / cm and at growth
rates down to 3pm/ min. Highly oriented Spherulites are obtained depending upon
the temperature gradients, degrees of superheat and pushing velocities employed. Y.
Huang [14] also conducted some experiments for polybutene—l samples. In the low
temperature region, the lamellar crystals grow in radial directions. However, once the
spherulites are large enough and contact each other, they would grow only into the
high temperature region caused by steric hindrance in the direction perpendicular to
the thermal gradient. Thus the application of a thermal gradient orients the lamellae

parallel to the thermal gradient within the Spherulite.

4.3 Previous work done in numerical modeling

Some of the earliest work done in modeling was by Tabar et al.[36] who Simulated
the two-dimensional Spherulitic growth and the effect of impingement on their small-
angle light scattering patterns. Galeski and Piorkowska [37] simulated the growth of
Spherulites for thermal, athermal and combined nucleation mechanisms in both two
and three dimensions. These simulations made it possible to examine the influence of
the nucleation mechanism on morphological details and verify the results of statistical
analyses of spherulite patterns. Some of the recent work in modeling microstructural
evolution was done by Billion et al.[38], Krause et al.[40] and Mehl and Rebenfeld [39].
These authors assumed that in a uniform thermal field the spherulites grow radially

outward from a spherical nucleus and used cellular discretization to track the growth

39

of a large pOpulation of spherulites and their impingement. Using this approach they
could predict the evolution of the microstructure. Billion et al.[38] correctly accounted
for the microstructural evolution of a thin film of neat polymer and the impingement of
Spherulites with the t0p and bottom surfaces. Similarly Krause et al.[40] modeled the
microstructure evolution of a polymer composite assuming that there are no obstacles
within the melt and that the entire domain is a neat polymer. The obstacles were
accounted for by deleting any volume of the spherulite that lies within an obstacle.
This limits their model as it cannot account for the “shadow” behind the obstacle.
Mehl and Rebenfeld [39] studied the solidification of polymer composites and modeled
the impingement between spherulites and inclusions by stopping growth for all radial
directions that intersect the obstacle. This approach is again applicable only when
the obstacles are much larger than the spherulites or when the obstacles traverse the
entire solidification domain. When these conditions are not satisfied, some errors
are introduced when it is compared to real systems. The most recent work done
was by Charbon et al. [44] who developed a numerical model for the solidification
of semicrystalline thermoplastic polymers. Their model is based on a front-tracking
method which accurately predicted the Shape of a spherulite but their approach was
computer intensive.

Since a number of these methods cannot account for situations where there are
small obstacles in the domain or where the growth velocity is a function of both
position and time, or the schemes used were computer intensive, it was decided to

use a model based on a standard front-tracking method that is efficient and accurate.

4.4 Basis of the Simulation

The Simulation model presented here for the growth of spherulites is based on solving

the two-dimensional Eikonal equation. The Eikonal equation relates the gradient of

40
the travel time to the velocity structure and allows Simulation of the propagation
of a two-dimensional wavefront for a given growth Speed. The eikonal equation is a
Special case of the Level Set equation which forms the basis of Level Set Methods.
They are numerical techniques that offer remarkably powerful tools for understand-
ing, analysing, and computing interface motion in a host of settings. Propagating
interfaces occur in a wide variety of settings. As physical entities, they include ocean
waves, burning flames, and material boundaries. These interfaces can develop sharp
corners, break apart, and merge together. The techniques have a wide range of
applications, including problems in fluid mechanics, combustion, manufacturing of
computer chips, computer animation, image processing, structure of snowflakes, and
the shape of snow bubbles. Rather than following the interface itself, the level set
approach takes the original curve and builds it into a surface. This surface is called

the Level Set function.

Let t(x,y) be the time at which a curve crosses the point (x,y). The surface t(x,y)

then satisfies the Eikonal equation

IVtIF =1, (4.1)

where F is the Speed of the front.

This equation indicates that the gradient of the time of arrival (t(x,y)) is inversely
proportional to the speed of the front. It is a so-called Hamilton-Jacobi equation when
the speed function, F, depends only on the position. Figure 4.3 Shows the outward
propagation of an initial curve and the accompanying motion of the level set function,
05. On the left a circle is Shown at time t=0 and below it is shown the circle at a
later time. In the upper right corner the initial position of the level set function is
shown and below we Show the same function at a later time. This is referred to as

the Eulerian formulation because the underlying coordinate system remains fixed.

Z : (X,y,l.=0)
l

h

---‘--——_
-_\-

 

 

 

z = (mat)

__-_—-__‘___--‘-
.-b
h-
‘-
‘

 

 

 

 

Figure 4.3: Illustration of a propagating circle and the corresponding zero level set

42

When expressed in Cartesian coordinates, the Eikonal equation takes the form

(is)2 + (ty)2 = 8207,10, (42}

where s(x,y) is the two-dimensional “slowness” (inverse of velocity or l/F) and
t:t(x,y) is the travel time field. This equation is a non-linear partial differential
equation which can be solved using finite difference schemes and a planar wavefront

extrapolation.

An optimum numerical method needs to be constructed to solve this partial dif-
ferential equation. Several methods have been developed to calculate travel times
on a regular grid. These methods or schemes are linked to those from hyperbolic
conservation laws. From a numerical point of view, the governing partial differential
equation can be approximated through the construction of an appropriate numerical
flux function, time ‘t’ in our case. The goal in the construction of such flux functions
is to make sure that the conservation form of the equation is preserved, to make sure
that the entropy condition is satisfied, and to try to produce smooth solutions away
from the discontinuities. The majority of these schemes are based on the idea of
viscosity solutions, which underly the conservation laws and upwind finite-difference
schemes. Therefore it is essential to clearly understand the concept of the viscosity

solution of Eikonal equations.

A very important property of viscosity solutions is that these are unique, weak
solutions, stable with respect to perturbations in the equation. Since roundoff and
discretization errors always perturb a numerical solution, the viscosity solutions are a
natural choice of numeric methods. However, the class of weak solutions is too big and
the initial data for these problems does not determine weak solutions uniquely. One
way to enforce uniqueness of the solution is to add viscosity to the Eikonal equation,

where viscosity is described by a dissipative term, which prevents the formation of

43
gradient discontinuities. Because smooth solutions are uniquely determined by the
appropriate initial data, the addition of viscosity fixes the uniqueness problem. The
notion of viscosity solutions was developed first for conservation laws, whose weak
solutions are made unique by a so-called entropy principle. The entropy solutions
of conservation laws are referred to as the viscosity solutions. Conservation laws are
commonly solved by upwind finite differencing, using either a backward or forward
finite—difference approximation to the derivative operator. Upwind finite-difference
methods are stable due to several reasons. Firstly, they add a minimal amount of
numerical viscosity to the equations, thereby tending to damp high-frequency error.
Secondly, there is no need to Specify Special numeric boundary conditions or rows of

fictitious points for these upstream schemes.

In the present work, two numerical schemes were used to solve equation 4.2. The
first one was based on the model developed by Vidale [45] to compute travel times
of seismic waves, while the second approach is based on recent work done by Sethian

[46].

4.5 Vidale’s scheme

The first scheme used in this work was based on the method devised by Vidale [45]
which computed travel times of the first arriving seismic waves on a two-dimensional
numerical grid by using a finite-difference extrapolation from point to point. The
method uses a velocity field that is sampled at discrete points in two-dimensional
space, with equal horizontal and vertical Spacing. Thus, the nucleus of the Spherulite
is assumed to be at the grid point A in Figure 4.4. The timing process is initiated by

assigning point A the travel time of zero. The four points adjacent to the point A,

44

 

   

A B1

 

 

 

 

I0 .93 .0
O

B4 C4

 

 

 

Figure 4.4: A schematic depiction of a source point A (nucleus) and its neighbors in
the ring

labeled as 81 through B4 are given the travel times
h
t, = 5(831-1" 8,4), (4.3)

where h is the mesh Spacing, 8A is the slowness at the point A, and 33,- are the
slowness at the grid points B,-. The formulae that extrapolate the travel times from
the three corners of a square to the fourth are derived below. Travel times between
the points A(t0 ), Bl(t1), and 82(t2) and the origin are sought. The first terms in the

equation are an average of two finite difference equations which are written as

 

 

0t _ to — t1

B—a: — h , (4.4)
and

9t- : t2 7 ’3. (4.5)

81: h

4.5

An average of these two equations gives

at 1
‘a—T Z Q—h-(to +t2 — t1 — t3).

Similarly, the second term is an average of the following two equations

 

 

:3:_t0—t2

8y_ h ’
and

gt_t1—t3

031— h '

The average of these two equations is

fl-ifl +t t—t)
8y_2ho 1 2 3-

Substituting equations (4.6) and (4.9) into equation (4.2) gives,

 

t3 = to + \/2(hs)2 — (I2 — t1)?

(4.0)

(4.9)

(4.10)

This equation gives the travel time to point C1 using the travel times from the source

to points A, B], and Bg. Similarly the arrival times are found for all the other corners.

Once all the four corners are found, the times are calculated for the entire ring and the

solution can progress to rings of increasing radius. The inductive scheme for adding

a ring of travel times to those already calculated is depicted in Figure 4.5. The

solution proceeds on the four Sides sequentially and is followed by the four corners.

The solution could start at any side arbitrarily, however one must take care that

the sides are traced sequentially. The critical point in these calculations is that the

travel times on any Side cannot be calculated in order. The points that are at a

46

 

0......
000000.
000000.
OOOOOO.
000000.
000000.
0......

O.

O.

 
 
 

 

 

 

 

Figure 4.5: A 2-D grid showing the numerical calculation of travel times

relative maximum and relative minimum need to be identified and the calculations
can proceed once those points are identified. A relative minimum is assumed if there
is a relative minimum in the time for the adjacent point in the adjacent row that has
already been solved in the previous ring. The minimum point on any edge is timed

using the plane-wave finite difference equation

 

t3 = to + \/(hs)2 — 0.250s2 — t1)2, (4.11)

where t3 is the time to be found, to is the relative minimum time in the inside row, and
t1 and t2 are the times on either side of the point whose time is to. Once the relative
minimum point is found, the solution progresses along the row finding the time for
each point until the relative maximum is encountered. Upon completion of the left
to right sweep, the row is swept in the reverse direction and the remaining untimed
points are solved in order from the relative minima to the relative maxima as shown

in Figure 4.6. This method is applied iteratively until the entire two-dimensional grid

47

 

Solvepoinlsatrehfivcmininn

Solvcbrkhtl’tunminimbmth

-———e- —.

9999900000..

Solvabldtfiunminimabmth

 

‘OOOOOOCDOOO

 

 

 

Figure 4.6: An illustration of the sequence of solution of one edge of the ring

is filled with travel times.

The principal drawback of this method is that the points on the computational
front must be ordered, with calculations starting at the so-called relative minimum
points. This ordering step makes the method computationally expensive. The choice
of the computational front is another limitation. Although the solution evolves in
time, it is computed by rectangular fronts that expand in Space. Thus if the time
field to be computed does not have an outward-pointing gradient at each front, i.e. if
the time gradient parallels the computational front, the square root in the governing
equations (4.10) and (4.11) becomes imaginary and the calculation must stop. This
situation occurs when the growth of Spherulites occurs in the presence of spatially
varying velocity. It fails when the ratio of these two ‘3’ values exceeds a certain limit.
The travel times in the ring just outside the obstacle boundary go complex because

of the Significant changes in the time gradients in and out of the obstacle.

4.6 Sethian’s approach

This second approach discussed here is based on the solution of the Eikonal equation

as discussed by Sethian [46]. This approach is schematically depicted in the flowchart

 

 

START
INITIALIZATION

 

 

 

 

 

END
INITIALIZATION

48

 

 

 

STOP ITERATIONS
PRINT THE FINAL
MATRIX.

 

 

 

 

 

 

BEGIN ITERATIONS
ON ALL THE NODES
IN THE GRID

 

 

 

 

CHECK FOR
CONVERGENCE

 

 

 

 

 

 

 

 

DEFINE FOUR
PARAMETERS OF
THE EQUATION

 

 

CONTINUE

 

 

ITERATIONS

 

 

 

 

 

 

SOLVE FOR THE EQN
AT EACH NODE.

REAL

 

 

 

 

 

CHECK FOR ALL
IMAGINARY VALUES

IMAG.

 

INTERPOLATE ‘T’
VALUES FOR
TAGGED POINTS

 

 

 

 

 

 

 

 

TAG ALL
THOSE POINTS

 

 

 

Figure 4.7: Flowchart for Sethian’s method

as Shown in Figure 4.7. A square grid of Size N x N was used to model the physical

Space. In this square box, the solution of the following equation is sought

mazI:(D,-;$t, O)2 + min(D$ft, 0)2 +

maxwgyt, 0)‘2 + mm(D:,iyt, 0)2 = F52, (4.12)
— __ 1: _tr— ,‘ +' _ tr .l_ti[ -y _ Iii—tw— +y __ in] —tI,l'
Where, Dijxt -- AI 1 , Dijxt — +1A1! , Dij t — Ay 1, Dij t — +Aly .

F is the speed function that is positive. F does not depend on the orientation of

the front and is simply a function of the applied thermal gradient. This is a quadratic

49
equation which is solved at each grid point. It is solved at each node over the domain
until convergence. These iterations are carried out several times through the entire

domain and the algorithm consists of:

(1) Initialization:- The source of the Spherulite which is the nucleus is considered
to be at the center of the grid. So for a N by N grid, the source point would be at
N/2 + 1, which is assigned the time value of zero. All the remaining points in the
box are then assigned random values which are increased away from the source point
in all directions in equal time steps. So if the origin lies at point (a,b) and the time
value at point (a,b) is ‘t’ (zero in this case), points in the positive x-direction would
be incremented as t + n x (increment), where “increment” is some fixed value and
n is the number of points. Once the time values for all the points on the x—axis are
initialized, all other points in the square box are timed similarly. Thus a square grid
is obtained filled with all time values which are initial guesses and which represent

the growth of a circular entity whose origin lies at the center of the grid.

(2) Solution:— The smallest time value was selected as a solution of the quadratic
equation for every point in the grid. The smallest value was selected to satisfy the so-
called entropy condition [46]. This condition ensures that the given surface will sweep
through the mesh in a consistent manner and that once the solution is obtained, it
is no longer reversible. This procedure of finding the smallest time values is carried
out throughout the domain starting at the lowermost left hand corner of the box.
At each iteration, the grid is scanned point by point until all the time values are
updated. The iterations are carried out several times until a converged solution is
obtained. During an iteration, complex solutions are sometimes encountered and
these are treated one by one using either an increment over the neighboring value or
simply taking an average value of the two neighboring points. If any such points lie

on either the x or y-axis, they are incremented over the neighboring value.

50
4.7 Implementation of the scheme

AS discussed earlier in section 4.4, this scheme is linked to hyperbolic conservation
laws. The governing equation (4.12) is simply an approximation of a hyperbolic con-
servation law through the construction of appropriate flux functions. For illustration,

the following hyperbolic conservation law is considered
II, + [C(u)], = 0. (4.13)

From the numerical point of view, the equation can be approximated through the

construction of appropriate numerical fluxes, g, such that

“In“ _ “I" _ .(J(ll'1nill‘i+ln) — g(ui—lnauin)
T _ _ A2: , (4.14)

 

where g(u,u) = C(u). One of the most straightforward approximate numerical fluxes

is the Engquist-Osher scheme[47], which is given by

”2 . dG
9E0(u17u2) = C(u1)+/ mm(

For G(u) = u2, the compact representation of this flux function is
ggo(u1, U2) 2 (7710.1‘('u1,0)2 + mz'n(u2, 0)2). (4.16)

This flux function is the core technique for approximating the level set equation

. . . . . I,~—I- -
under conSideration. In this problem, for the one-dimenSIOnal case, ul :2 #1,

t

I —t‘ . . . . . . 0
U2 = 4153—1. This scheme works mainly because it satisfies the condition of conSIS-

tency and also because each grid point in the domain is updated using the smaller

neighboring values. This provides a consistent way of sweeping through the mesh

51
and constructing the solution surface ‘t’. The same technique can be extended to the
two or three dimensional case. This scheme though very versatile and more generic
has certain limitations. It is found that this scheme is extremely sensitive to initial
conditions. The validity of the results depend on the way the process of initialization
has been done. The Simulations were carried out using a variety of initial conditions.
Most of the initial conditions tried resulted in a large number of complex numbers
and the time values could not be interpolated. It was found that the initialization
procedure as discussed in section 4.6 is the most efficient way and the results obtained
were very reasonable. Therefore extreme care needs to be taken to assign the initial

time values.

4.8 Results and Discussion

The simulations presented below are carried out using Sethian’s approach. Though
Vidale’s scheme gives correct results for single and multiple spherulites, it fails to
account for the growth in presence of fibers. The travel times go complex because of
changes in the time gradient when the spherulite encounters an obstacle as discussed
earlier. Hence the simulations were carried out using Sethian’s method as it is versatile
and has no limitations. The method is relatively simple and can be used for very large
grid sizes. Also, if there are any complex numbers encountered, they can be treated

easily.

4.8.1 Spherulitic Growth in a Neat Polymer

Figure 4.8 Shows the shape of a spherulite predicted under isothermal conditions
using the above scheme. The shape of a Spherulite is a circle which has been observed
by many experiments [1]. This simulation has been carried out under a uniform

melt temperature condition which indicates a constant growth velocity. Thus all the

52

 

 

 

’ ................ x
I’ _’ ‘-.
............ .-
l‘ ‘-
/' \-
o \.
l./ \.
. ’ ..... ‘- \ -
I . I ' \ .\
I .’ ‘ '
I ’I \ \.\
'I I- ..... \\ -\ '. -I
I - . '.
. I \ \ .
I ' . 1.
. I \ \ .
I . \ &
. I ‘ \ ' '.
I I 0 . \ o
‘ I ' i
. , ‘
[ I . 3 .
I I I "
\ It ,‘ I [
‘ \ ,- ! .=
'\ .\ ' .l f
'\ \ .’ I ,.-'
.\ \ I ‘I I :-
.\ \ ~- -. ” .I -
Q ' I
\ ’ .
\, ~. , /
\ ‘_ ‘ ’ ’ I
\. ' ‘ ' ' _/. —I
\.\ I
.‘.‘ ’I" o
---------- I’ u.
.......... "
I l I """" L. I 4.1 ----- I I " I l

 

10 20 30 40 50 60 70 80 90100

x (microns)

”final.dat" --

70 --
60
50 .....

100 40 ~-
30 -~

90 20
10 -

80

70

60

50 y (microns)

40

30

20

10

Figure 4.8: Front position of a spherulite at different times while growing in an
uniform temperature field (Grid size = 100x100 microns and the times are in mi-

crons/sec)

 

"iinal.dat" —-
so
._.--" .--“. ,-"‘. .‘*-_. 1 100 30 ---
L'r‘ ‘. .~““ 20 -.
..................... ........................ ._ ".f 90 1o
I ————— ~"\_ .’ I """"" ‘. .1 80
I, \,\ / \\
’ \1
I , ’ ' ‘. ..... ‘ \\ ‘1 70
I /. \\. ’.r _\ \.
I‘ \‘ .’. .\ \
I. ', .. \_l .' ‘. \. _‘
’ .' ‘. \, .' .. ‘ , .
i ~. : i -. : 3 2 - 50 y (microns)
" ‘. - .' I .""\ ‘. . .' ' .i I
\‘\ ~.___.- ./ .\ ~.._-.- '1’ [I - 40
- / .‘ ’ '
\ I . . I
\. -‘-‘ ’." ‘.~. l’,’ ’-
\'\ “'-' .x -m’ .’. 7‘ 30
'\ / 'x / ..‘
'\ ~ I I ‘\.\ ’ I 20
....... I -‘_‘—.’.’ . .- q
‘<. ................. '0 h ...................... --' "..'1 10
"In- I 1 1"“) l 1 1'“ 1 ,

 

 

 

102030405060708090100

x (microns)

Figure 4.9: Shape of two spherulites at different times illustrating impingement and
subsequent growth (Grid size = 100x 100 microns and the times are in microns/sec)

spherulites grow with the same velocity. Since the simulations give satisfactory results

for single spherulites, we now extend the calculations to two or more spherulites.

Figure 4.9 indicates the growth of two spherulites whose nuclei lie at a certain
distance apart in the domain and the solidification progresses at a constant growth
velocity. As time progresses, the growth fronts from neighboring spherulites impinge
on each other. Figure 4.10 is another such example where the nuclei lie at an off-

centered location in the domain.

'final.dat" --

30 ...

25 -----

, ,. .. fl ..... _ ~. 100 20 -..

 

V ' .-""'N‘"...' /'—"-.. ' ._.
,v -- ,~ ~. .-.-‘ ‘- 1 -.
, ..‘ ” \ I .‘ '. s
. ” , ,s '\_l . ‘- \
, ' ‘ -- ~ A 90 10
I. "‘ _/ ._. \_ I x. \ ._.
/- I. ‘. I ,' F'. \ . “ - '.
I " ’.".~ I ' I- u. .' s .

.
o \ ‘
r .' '.. | .n . ' . I '5 ‘ 1 . !
o a . \ o , | ‘ .
| , ' I n “ '.. ' , ' ‘ I c I I
, .
|. . - " . '. .‘ .‘ x a ,.' . I' I J 70
\ \ ‘. ‘~ . \ ‘- -. —’ t. -, x . II _.'
. . . ' ‘ \ .. ,v I '
. \ \ . . . ' .I s, ’. . , I. ’
-. - . ~ ., - \ ~. . \\ ’. _ ..
‘. \ s , I x. . ‘ i . I "- _ w
'. . . ‘,...‘ \ . ,... .- . ,-
: -/ -, ' ‘. ‘ I " ~ ‘. \ ,I" ." 0‘

i 5 a" z" ‘2 '3 {' z" '.= .3 3 i 4 50 y(microns)

0".
u
-.

I v ’ . . ' .'. \
u " /. \ g 'I .\ '.. J
- \ ‘--.- . ‘~ ..' ' .. . ‘- 4”
.\ ' .- . 'I .l " ‘-\ \. '..
O. . ~‘ \- .0.

 

.' I I I , .’
-.' . I 4 ‘ ‘ o 0" .. \. .’ ‘ a“\. ‘. \u u.
. I I , , \ a u v v ~ \ \ '.
' . . , .’ ‘~ . o ' I ‘ . . -
I I ' g, ‘. ‘ u. v-..‘ “ ,' o'.‘\ X \ \ 1 30
Z ' . ' u n ' v r “ . o I ‘ . ‘ '
'. ‘. \ '. ~. ,’ .‘ - i . ' ~ : ‘. ~ . !
-' . ‘ ._. ' I ‘ . a o ‘ I a
\_ \. -. .. — 1' “ ‘ .' I .’- ‘. ‘ " , _I ,_ ._.
. -_. \ \‘ 1' \ ‘-- I \ *.. .I ' ;
‘ \~ \ ‘ ‘ I, ‘. .’ .’ .’
‘ \. ‘- ‘-_." ‘___,r I . _ 1o
, \‘ .l' \ .A- . .‘I .
L . ".' '\ o’ ‘ o/ 0' I
‘ ......... ‘ ‘‘‘‘‘ ’
‘x l m I """ I' ha. I 4-1 .1 L ---- I

 

 

102030405060708090100

x (microns)

Figure 4.10: Impingement of multiple spherulites with randomly located nuclei shown
at different instants (Grid size = 100x 100 microns and the times are in microns/sec)

55
4.8.2 Spherulitic Growth in the Presence of Obstacles

The growth of spherulites is also studied in the presence of obstacles which are often
present in composites in the form of fibers. The addition of fibers to a given crystal-
lizing system increases the complexity of the system by adding fiber surfaces which
can both constrain growth by an impingement mechanism and enhance crystalliza-
tion by adding nucleation sites. In a first example, a single spherulite growing in the
presence of a square inclusion is shown in Figure 4.11. This figure shows that there is
an excellent agreement between these results and those obtained by Schulze and Bier-
mann [42] in their experiments. The simulations exactly reproduce the experimental
shape of the spherulite for all the time steps. The growth morphology maintains a
front normal to the surface of the inclusion as observed experimentally. Figure 4.12
is another case where the obstacles are assumed to be rectangles and randomly dis-
tributed within the solidification domain. The modeling is based on the fact that
when the spherulite encounters an obstacle, the lamellae that are in contact with the
fiber stop growing while the rest of the lamellae continue to grow at the normal speed.
Figure 4.12 shows that the simulations give satisfactory results for the presence of
multiple obstacles in the domain. For the implementation of the numerical scheme,
when an obstacle is encountered, the nodes or points that lie in the region marked as
an obstacle, have their local velocity reduced considerably as compared to the growth
velocity elsewhere in the solidification domain (typically 1% of the growth velocity).
The growth front eventually wraps around the obstacle and a microstructure which

maintains a front normal to the surface of the inclusion is obtained.

56

“final.dat" --

7o --

60

so

.. 100 40
........... 30
v ’ ~ 90 20
.1" ‘.-.-.‘. \\ 10 "-

 

I a ‘ ‘ . \ ‘.
’l' ’.,- s.\ .\ _‘ 80
.’ /' .‘. \

.’ I' . ----- ~. ‘. '\

I ‘ x ‘- . \ ‘ ‘.
, .I x -‘ -\ \. 1 7O

. / I ~ , \
I _ ' \ '\

I!- f, .v "‘ a. a“ K.“ K. .‘ fi 60
l ; ; l :y s 50 y (microns)
k, , ,5 - 4o

' .\\. . . ,. x ,- " /. .’. 7' 30

 

 

 

 

\ p I. I

I. ‘. ,

\' \‘ . . ’ I ‘ I "/ _.'.. d 20
) \ ‘ - ~ . - - ' .l l:
'\ .’

"- x" .4 10

\'\ .......... ’ _.".. 0'

'

 

 

‘-
a
o
n'
.
o

 

102030405060708090100

x (microns)

Figure 4.11: Growth of a spherulite in the presence of a rectangular obstacle shown
at different times (Grid size = 100x100 microns and the times are in microns/sec)

57

"final.dat" --
7O --

 

",r f ‘ 100 40 ---
1”." \N “X.“ - 90 20 "-2
. i h -------- q 30
i .I i "N. inf-— 70
.I' ,' .. '~ -'\. - 60

‘ ,' .‘ h ‘.

 

 

 

 

 

i . i j E- 50y(microns)

'\ '\. ' I." ." ,I' ‘ 40

 

 

 

 

 

\ ‘‘‘‘‘ .’ I
f \'\ " ‘ .’ I. 30
'._ . . " u, / _’l #
\ \ -. __ _,' I. .I
T. " _I. .l
\. ‘\ i I ' I I, d 20
\. '~ ....... ‘,-
\'\ ."’
.\'\ ......... "I! on"! .'- 10
.-" I.

 

1O 20 30 40 50 60 70 80 90100

x (microns)

Figure 4.12: Spherulitic growth in presence of multiple obstacles shown at different
times (Grid size = 100x100 microns and the times are in microns/sec)

58

4.8.3 Spherulitic Growth in the Presence of a Temperature

Gradient

Finally, a spherulite that grows in a temperature gradient is considered. This is
done to simulate, for example, directional solidification or growth near the wall of
a mold where high gradients are present. A linear temperature gradient (of the
order of 104°K/ m) was applied and the shape of the spherulite obtained from the
simulations closely matches the morphological features appearing in the polarization
photomicrographs from the work done by Lovinger et al.[43]. The nucleation rate
is sufficiently suppressed and the fibrils of the growing spherulite curve towards the
direction of motion of the crystallizing front. Figure 4.13 shows this case of the
curving of spherulite fibrils from left to right. This shape resembles the teardrop
shape obtained in the experiments performed earlier by Lovinger et al.[43], which is

the result of the different growth velocities of the solid and liquid phases.

4.9 Conclusions

Two different methods have been discussed in this work which can simulate the growth
of spherulites from the melt as a function of the conditions present during the actual
solidification process. Both techniques are based on solving the Eikonal equation
to track the growth front. Though the first scheme used for modeling is accurate
for simulating growth of single or multiple spherulites, it fails to model the growth
when there are differences in the local velocities of points due to the appearance of
complex numbers during the calculations. The second method proposed by Sethian
is more versatile and can predict the time evolution of the shape of a spherulite. The
technique is robust and is based on finite difference computation which can be used for

a large grid size as well as for growth of multiple spherulites. It can also address the

59

"final.dat" --

50 --
30 .....
, ,. .. 100 20 ---

 

 

 

/ . .
’1’ '5‘ ......................... ff‘. 10 ....
x ------- - 90
, ........ - so
,I' ‘x. ".
I" \x
f, \'\, T 70
/ . - - a \
/ . / ' \ . * 1'
/' ,-’ ‘. \ '-.' 60
I. ,' \. "
' . \ . ': .
i I. .. y 4 50 y (microns)
\ \ I .' .5
. - ’ I :-
‘.\ ‘.\‘ . ’.’ _I ::_ 40
. 1. .\ x _ _ . / I ,-
1. .\ .I
\‘ .\. ./'/ ‘ 30
1 -, \. oz
‘s ..'«‘ \ ° \ ,,,, ' I 3'.-
, "‘ _ 20
\ \\ .'..- "."
\ ‘ .-'
\ .‘ v""o. ,.""'
\\ ‘ . ......................... i .f‘ 10
\\ s‘ '. '
B 1 '1 1 1 1 1 n 1

 

10 20 3O 40 50 60 70 80 90100

x (microns)

Figure 4.13: Spherulitic growth in the presence of a temperature gradient shown at
different times (Grid size = 100x 100 microns, times are in microns/sec and dT/dx
= 104°K/m )

60
problem of the growth of spherulites in the presence of fibers and in the presence of an
applied thermal gradient. This was demonstrated with a variety of two—dimensional
examples. Both these methods take approximately the same computational times on
similar platforms. The simulations were performed on Sun UltraSparc with a CPU
speed of 167 MHz and 64 MB Ram. The operating system was SunOS Release 5.6
(UNIX System V Release 4.0). The domain size for the all the results presented
was 100x 100 microns. The approximate time for each simulation performed was 300

seconds.

Chapter 5

Summary and Conclusions

The solidification of semicrystalline polymers is studied in this work using numerical
and analytical models. More specifically, the work presented is directed towards the
possibility of using directional solidification to obtain polymer/ polymer composites
and to model the evolution of the microstructure.

A brief review of the literature pertaining to the types of polymers, and crys-
tallization in semicrystalline polymers was first presented. A numerical scheme to
describe the evolution of spherulites under various conditions was then presented
Two different methods were employed to simulate the growth of spherulites, Vidale’s
scheme and Sethian’s approach. Both these techniques were based on solving the
Eikonal equation. The method proposed by Sethian is more versatile and can predict
the time evolution of single and multiple spherulites, growth of spherulites in the
presence of fibers and also in the presence of an applied thermal gradient. Several
simulations were carried out to demonstrate these methods. The cases of isother-
mal crystallization, crystallization in a temperature gradient, and crystallization in
the presence of obstacles are considered and various examples are solved for two-
dimensional processes. Crystallization of a homogeneous polymer blend composed of

a non-crystallizable polymer and a crystallizable polymer was then considered. An-

61

62
alytical work was done directed towards studying the possibility of forming in-situ
polymer/ polymer composites using directional solidification. The vital stages of the
non-isothermal processes of microstructure formation such as nucleation and growth
were discussed. The thermal conditions that favor steady growth and which suppress
nucleation ahead of the growing interface were identified. The stability issues of the
growing interface were then addressed. It was shown that for the given problem un-
der consideration, the interface always became unstable, which gave rise to crystalline
fingers in the direction of the thermal gradient. The work presented here is hoped
to serve as the basis for future engineering designs of materials and processes related
to this new type of composite. Shrinkage and drawability of these new materials are
expected to be very good, with improved mechanical properties compared to the pure
amorphous polymer. It is hoped that this work will enhance our understanding of the
crystallization process in polymers. Also, the numerical model presented is a versatile
model describing the shape of polymer spherulites formed in a temperature gradient.

The computational times are very short as compared to other models developed.

Chapter 6

Suggestions for Future Work

In this work, an efficient two-dimensional numerical model has been developed that
satisfactorily predicts the microstructure and crystallization kinetics of polymer com-
posites. However it would be more realistic to develop a model that can perform
three-dimensional simulations. The present model considers instantaneous nucleation
in an isothermal field. Future work could involve the study of the influence of the
mode of nucleation and the functionality of the Spherulitic growth rate. The numerical
model presented here simulates the growth of spherulites in the presence of obstacles,
but it does not consider the types and size of fibers. A close examination could be
made of the effect of fiber size and shape on the crystallization process in reinforced
systems for both two and three dimensional cases. Once the three-dimensional model
is accurately deve10ped it can be coupled with a finite element code developed. A
model based on an enthalpy formulation of the energy equation can simulate the so—
lidification process occuring in complex three-dimensional parts when combined with
these simulations. It will thus be possible to compare the observed morphology of a
complex part with the computed results.

It was also found in the course of performing the simulations that the initial

conditions have some influence on the rate of convergence of the numerical scheme.

63

64
A combination of both Vidale’s and Sethian’s schemes might provide a methodology

that will improve the convergence of the algorithm.

APPENDICES

Appendix A

C code based on Sethian’s

approach

'/.'/.'/.'/.'/.°/.'/.'/.'/.'/.°/.'/.'/.'/.'/.°/.°/.'/.°/.'/.'/.'/.'/.'/.'/.'/.'/.'/.°/.'/.°/.'/.'/.'/.'/.'/.'/.'/.'/.'/.7.7.'/.'/.7.'/.%'/.°/.'/.'/.'/.'/.'/.'/.'/.7.7.7.'/.°/.°/.'/.7.7.7.
Z This is the main simulation program based on Sethian’s approach%
'/.'/.'/.'/.'/.°/.'/.'/.'/.'/.'/.'/.°/.°/.'/.'/.'/.°/.'/.'/.°/.'/.'/.'/.°/.'/.'/.'/.'/.'/.°/.°/.°/.'/.°/.°/.'/.'/.'/.'/.'/.°/.'/.'/.°/.'/.'/.'/.'/.'/.'/.'/.'/.'/.'/.'/.'/.'/.°/.'/.'/.'/.'/.'/.°/.'/.
#include<stdio.h>

#include<math.h>

double TESOOJESOO], TNEW[300][300];

double F[300][300];

double slowF[300][3001;

int N = 200;

double dX, dY;

float initVal = -0.05;

float finVal;

float orgVal = 0.0;

float tol =0.0000010;

float sum;

float incr;

65

66

int obsLen, obsBrt, obsOrgX, obsOrgY;
struct _obj
{
int obsLen;
int obsBrt;
int obsDrgX;
int obsOrgY;
}
object[100];
double min(doub1e x, double y)
{
if(x<y) return x;
else return y;
}
double max(doub1e x, double y)
{
if(x>y) return x;

else return y;

main()

int i, j, a, b, k, counter;

int flagl, flag2, f1ag3, f1ag4;
float fact1, fact2, fact3, fact4;
double A, B, C;

/* double xdim=100; */

double rootVal;

double p,q;

67
int iter=0;
int n_objects;
int done=0;

int centX[5], centY[5], cents, N0_CENTS;

incr = 10.0;
dX = dY = 1.0/N;
for(i=0; i<N+2; i++)
for(j=0; j<N+2; j++)
{
F[i][j] = .01;
slowF[i][j] = 0.000001;/*1/F[i][j];*/
}
printf("Give no. of objects\n");
scanf("%d",&n_objects);
for(i=0; i<n_objects; i++)
{
printf("Give obs. length\n");

scanf("%d",&object[i].obsLen);

printf("Give obs. breadth\n");

scanf("%d",&object[i].obsBrt);

printf("Give obs. xorigin\n");

scanf("%d",&object[i].obs0rgX);

printf("Give obs. yorigin\n");

scanf("%d",&object[i].obsOrgY);

68

printf("Data inputed:\n");
for(i=0; i<n_objects; i++)
{
printf("Zd ",objectEi].obsLen);
printf("Zd ",objectfi].obsBrt);
printf("xd ",object[i].ob80rgX);

printf("%d\n",object[i].obsOrgY);

printf("print the centers, print -1 -1 to stOp : \n");
ND_CENTS=0;
for(i=0;i<5;i++)
{
scanf("%d Zd", &centX[i], &centY[i]);
if (centX [i] ==-1) break;
N0_CENTS++;
}
for(cents=0;cents<N0_CENTS;cents++)
{
a = centX[cents];

b = centYEcents];

T[a][b] = orgVal;

for(k=b+1; k<=N+1; k++)

T[a][k] = T[a][k-1]+incr;

for(k=b-1; k>=0; k--)

T[a][k] = T[a][k+1]+incr;

for(i=b; i<=N+1; i++>

for(j=a+1; j<=N+1; j++)
{
T[j][i] = T[j-1][i] + incr;
}
for(j=a-1; j>=0; j--)
{

T[j][i] = T[j+1][i] + incr;

for(i=b; i>=0; i--)

for(j=a+1; j<=N+1; j++>
{
T[j][i] = T[j-1][i] + incr;
}
for(j=a-1; j>=0; j--)
{

69

70

T[j][i] = T[j+1][i] + incr;

}
}
}
while(1)
{
for(j=1; j<N+1; j++)
{

for(i=1; i<N+1; i++)
{

if( (i==a) && (j==b) ) continue;

factl = (T[i][j] - T[i-1][j])/dx;
fact2 = (T[i+1][j] - T[i][j])/dx;
fact3 = (T[i][j] - T[i][j-1])/dY;
fact4 = (T[i][j+1] - T[i1[j])/dY;

flag1=f1ag2=fla33=f1ag4=1;

if(fact1<=0.0) {factl

0.0;f1agi=0;}

if(fact2>=0.0) {fact2 0.0;f1ag2=0;}

if(fact3<=0.0) {fact3

0.0;f1ag3=0;}

if(fact4>=0.0) {fact4

0.0;f1ag4=0;}

b
ll

(f1ag1+f1ag2)*dY*dY + (flag3+flag4)*dX*dX;

on
ll

-2*( dYtdY*(T[i-1][j]*f1ag1+T[i+1][j]*flag2) +

71

dXth*(T[i][j-1]*flag3+T[i][j+1]*flag4) );

/* if( (i>obs0rgX) && (j>obsDrgY) ) */
for(counter=0;counter<n_objects;counter++) {
obsLen = objectfcounter].obsLen;

obsBrt = object[counter].obsBrt;

obsOrgx object[counter].obsOrgX;

obsDrgY object[counter].obsOrgY;

if( (i>obsOrgX) && (i<obsOrgX+obsLen) && (j>obsOrgY)
&& (j<obsOrgY+obsBrt) ) {
done=1;

break;

if(done)

C)
II

dYtdY*(T[i-1][j]*T[i-1][j]*f1ag1+T[i+1][j]*T[i+1][j]*flag2)
+ dXth*(T[i][j-1]*T[i][j-1]*f1ag3+T[i][j+1]*T[i][j+1]*f1ag4)
- dXth*dY*dY/(slowF[i][j]*slowF[i][j]);

else

0
ll

dY*dY*(T[i-1][j]*T[i—1][j]*f1agi+T[i+1][j]*T[i+1][j]*f1ag2)
+ dXth*(T[i][j-1]*T[i][j-1]*f1agS+T[i][j+1]*T[i][j+1]*f1ag4)

- dXth*dY*dY/(F[i][j]*F[i][j]);

done=0;
if ( (A == 0.0) II ((B*B-4.0*A*C) < 0.0) )

TNEw[i][j] = -3.00;

72
else
TNEW[i][j] = (-B+sqrt(B*B-4.0*A*C))/(2.0*A);

}

for(i=1; i<N+1; i++)

for(j=1; j<N+1; j++)
{

if(TNEW[i][j]!=-3.00) continue;

/*find quadrant*/

if( (i>a) && (j>b) > /*lst*/

p = max(TNEw[i-1][j]. TNEWIiJEj-ll):

q = min(TNEw[i][j+1], TNEW[i+1][j]);

if(p==-3.0)
printf("p==-3.0\n");
else if(q==-3.0)
TN£w[11[j] = p + incr;
else if( (i==N) ll (j==N) )
TNEw[i][j] = p+incr;
/*else if(p>=q)
TNEHfi][j] = p+incr;*/

else

73

TNEWEiJEj] = (p+q)/2.0;

else if( (i<a) && (j>b) ) /*2nd*/

p = max(rNEwEi+1][j]. TNEWEi][j-1]);

q = min(TNEH[i][j+1]. TNEWfi-llfjl);

if(p==-3.0)
printf("p==-3.0\n");
else if((TNEw[i+1][j]==-3.0) II (TNEW[i][j+1]==-3.0))
TNEW[i][j] = (TNEW[i][j-1] + TNEWEi-l][j])/2 0;
else if( (i==1) ll (j==N) )
TNEHfilfj] = p+incr;

/*else if(p>q)

TNEH[i][j] = p+incr;*/
else
TNBWEilfj] = (p+q)/2.0;
}
else if( (i<a) ea (j<b) ) /*3rd*/
{

p = max(TNEW[i+1][j], TNEW[i][j+1]);

q = min(TNEW[i-1][j]. TNEw[i][j-1]);

if(q==-3.0)
printf("p==-3.0\n");

else if(p==-3.0)

74
Tnsw[i][j] = q + incr;
else if((TNEW[i+1][j]==-3.0) II (TNEW[i][j+1]==-3.0))
TNEWfilfj] = q + incr;
else if(i==1)
TNEUEi][j]=TNEH[i-1][jJ-incr;
else if (j==1)
rnsw[i][j]=TNEw[i][j-1]+incr;

else if(p>q)

TNEW[i][j] = p+incr;
else
TNEW[i][j] = (p+q)/2.0;
}
else if( (i>a) && (j<b) ) /*4th*/
{
p = max(TNBw[i][j+1], TNEWEi-1][j]);
q = min(TNEW[i][j-1], TNEWEi+1][j]);
if(p==-3.0)

printf("p==-3.0\n");
else if(q==-3.0)

TNEW[i][j] = (TNBw[i][j-1] + TNEw[i-1][j])/2.o;
else if( (i==N) || (j==N) )

TNEWfi][j]

p+incr;
else if(p>q)

TNEH[i][j] = p+incr;
else

rnswtiltj]

(p+q)/2.0;

75

else if(j==b) /*x axis*/

if(i>a) /*+ve*/
TNEW[i][j] = TNEWEi-1][j]+incr;
else

TNEw[i][j] = TNEWEi+1][j]+incr;

else if(i==a) /*y axis*/

if(j>b) /*+ve*/

TNEVEi][j] = TNEW[i][j-1]+incr;

else

TNEWEilfj] = TNEWEi][j+1]+incr;

else

printf("UUTTA MY HANDS\n");

/* west */

if(i==1)

TNEW[i-1][j] = TNEwti+1][jJ;

76

/* east */

if(i==N)
{
TNEH[i+1][j] = TNEw[i—1][j];
}
/* n*/
if(j==N)
{
TNEWEi][j+1] = TNEW[i][j-1];
}
TNEWEa][b1=0.0;
}
}

/*compare T & TNEW*/

sum=0.0;

for(i=1; i<N+1; i++)
for(j=1; j<N+1; j++)

sum += (TEiJEjJ-TNEwti][j])*(T[i][jJ-TNsw[i][j]);

/*if NOT 0K, copy T to TNEW */

printf("DIF = %f\n", sum);

if(sum <= tol) break;

for(i=1; i<N+1; i++)

77
for(j=1; j<N+l; j++)

T[i][j] = Tnsw[i][j];

iter++;
if(iter>600) break;

} /*end of while*/

finale :
printf("The final Matrix\n");
for(i=1; i<N+1; i++)
{
for(j=1; j<N+1; j++)
printf("xz.2f ", TNEw[i][j]);
printf("\n");
}

printf("# of iters=%d\n", iter);

Appendix B

Matlab program for simulations

based on Vidale’s scheme

B.1 mex.m

ZZZZXZXZXXXZZXZZZZXZXXXZZXZZXXZXXZXZZZZZ1ZZZZZZZZZXZZZZZZZZZZX
% This is the main program which inturn calls the other files%
% chmi.m, middle.m, corner.m, Z
%%%%%%%%%%ZZZZZZZXZZXXZZZXZZZZZZZZZZXZXZZXZZZZZZZZZXZXXXZZZZXX
function [valuecorn]=corner(t0,t1,t2,xpo,ypo)
global ring sid post hsori calmid rm cm gr hsobs xbeg
xfin ybeg yfin
global obx oby xsizob ysizob inflec cond least
hs=hsori;

if(xpo>=obx+1)

if(xpo<=obx+xsizob-1)

if(ypo>=oby+1)

if(ypo<=oby+ysizob-1)

hs=hsobs;

78

79

end
end
end
end
ring(xpo,ypo)=t0+sqrt((2*(hs“2))-((t2-t1)“2));
prevrin=t0+sqrt((2*(hs‘2))-((t2-t1)“2));
valuecorn=t0+sqrt((2*(hS‘2))-((t2-t1)“2));
if(gr>1)
for chek=1z10
if(sid==1)
if(least(4*(gr-2)+1,chek)==xpo)
t0=ring(xpo,ypo-1);
t1=ring(xpo+1,ypo-1);
t2=ring(xpo-1,ypo-1);
ring(xpo,ypo)=t0+sqrt(((hs‘2))-.25*((t2-t1)‘2))3
valuecorn=t0+sqrt(((hs“2))-.25*((t2-t1)“2));
sidel=1;
xpo=xpo;
YP°=YP°3
Zpause
end
end
if(sid==2)
if(least(4*(gr-2)+2,chek)==ypo)
t0=ring(xpo+1,yp0);
t1=ring(xpo+1,ypo+1);
t2=ring(xpo+1,ypo-1);

ring(xpo,ypo)=t0+sqrt(((hs“2))-.25*((t2-t1)“2));

80
valuecorn=t0+sqrt(((hs“2))-.25*((t2-t1)“2));
side2=2;
xpo=xpo;
ypo=ypo;

Zpause
end
end
if(sid==3)
if(least(4*(gr-2)+3,chek)==xpo)
t0=ring(xpo,ypo+1);
t1=ring(xpo+1,ypo+1);
t2=ring(xpo-1,ypo+1);
ring(xpo,ypo)=t0+sqrt(((hs“2))-.25*((t2-t1)‘2));
valuecorn=t0+sqrt(((h8“2))-.25*((t2-t1)“2));
side3=3;
xpo=xpo;
YP°=YP°i
Zpause
end
end
if(sid==4)
if(least(4*(gr-2)+4,chek)==ypo)
t0=ring(xpo-1,ypo);
t1=ring(xpo-1,ypo+1);
t2=ring(xpo-1,ypo-1);
ring(xpo,ypo)=t0+sqrt(((hs‘2))-.25*((t2-t1)‘2));
valuecorn=t0+sqrt(((hS“2))-.25*((t2-t1)“2));

side4=4;

81
xpo=xpo;
YP°=YP°3
Zpause
end
end
end

end

B.2 chmi.m

function [valuecent]=chmi(gri)
global ring rm cm post sid inflec cond least
smal=1e37;
for side=1z4
smal=1e37;
for thu=rm-gri:rm+gri
if(side==1)

if(ring(thu,cm+gri)<smal)

sma1=ring(thu,cm+gri);

post(gri,side)=thu;
end

end

if(side==2)

if(ring(rm-gri,thu)<smal)

snual=ring<rm-gri,thu);

post (gri ,side)=thu;

end
end
if(side==3)

if(ring(thu,cm-gri)<smal)

sma1=ring(thu,cm-gri);

post(gri,side)=thu;
end

end

if(side==4)

if(ring(rm+gri,thu)<smal)

smal=ring(rm+gri,thu);

post(gri,side)=thu;
end

end

end Zfor thu

z-_----__-_-________-_--_-____-_; ___________________________

if(side==1)

remup=(rm-gri);
remdo=(rm+gri);
for thu=post(gri,1):-1:remup+1

if(ring(thu,cm+gri)>ring(thu-1,cm+gri))

83

ind=4*(gri-1)+1;
cond(ind)=cond(ind)+1;
inflec(ind,cond(ind))=thu;
least(ind,cond(ind))=thu—1;
sidelU=1;
grid=gri;
Zpause

end

end

for thu=post(gri,1):remdo-1

if(ring(thu,cm+gri)>ring(thu+1,cm+gri))

ind=4*(gri-1)+1;
cond(ind)=cond(ind)+1;
inflec(ind,cond(ind))=thu;
1east(ind,cond(ind))=thu+1;
side1D=1;

grid=gri;

Zpause

end
end
end choses side=1

Z ___________________________________________________________

if(side==2)

84

remlft=(cm-gri);
remrgt=(cm+gri);
for thu=post(gri,2):-1:remlft+1

if(ring(rm-gri,thu)>ring(rm-gri,thu-1))

thu=thu;

first=ring(cm-gri,thu);

second=ring(cm-gri,thu-l);
ind=4*(gri-1)+2;
cond(ind)=cond(ind)+1;
inflec(ind,cond(ind))=thu;
least(ind,cond(ind))=thu-1;
side2L=2;
grid=gri;

%pause

end
end
for thu=post(gri,2):remrgt-1

if(ring(rm-gri,thu)>ring(rm-gri,thu+1))

ind=4*(gri-1)+2;
cond(ind)=cond(ind)+1;
inflec(ind,cond(ind))=thu;
1east(ind,cond(ind))=thu+1;
aside2r=2;
E;rid=gri;

'Lpause

end
end

end choses side=2

% ___________________________________________________________

if(side==3)

remup=(rm-gri);
remdo=(rm+gri);
for thu=post(gri,3):-1:remup+1

if(ring(thu,cm-gri)>ring(thu-1,cm-gri))

ind=4*(gri-1)+3;
cond(ind)=cond(ind)+1;
inflec(ind,cond(ind))=thu;
least(ind,cond(ind))=thu-1;
side30=3;

grid=gri;

Zpause

end
end
for thu=post(gri,3):remdo-1

if(ring(thu,cm-gri)>ring(thu+1,cm-gri))

ind=4*(gri-1)+3;
cond(ind)=cond(ind)+1;

inflec(ind,cond(ind))=thu;

86
least(ind,cond(ind))=thu+1;
side3D=3;
grid=gri;

Zpause

end
end
end Xcloses side=3

Z ___________________________________________________________

if(side==4)

remlft=(cm-gri);
remrgt=(cm+gri);
for thu=post(gri,4):-1:remlft+1

if(ring(rm+gri,thu)>ring(rm+gri,thu-1))

ind=4*(gri-1)+4;
cond(ind)=cond(ind)+1;
inflec(ind,cond(ind))=thu;
least(ind,cond(ind))=thu-1;
side4L=4;

grid=gri;

Xpause

end
end

for thu=post(gri,4):remrgt-1

87

if(ring(rm+gri,thu)>ring(rm+gri,thu+1))

ind=4*(gri-1)+4;
cond(ind)=cond(ind)+1;
inflec(ind,cond(ind))=thu;
1east(ind,cond(ind))=thu+1;
side4R=4;

grid=gri;

Zpause

end
end
end choses side=4
Z ___________________________________________________________
end Zfor side
post;
Zpause

valuecent=10;%You know

B.3 middle.m

function [valuemid]=middle(t0,t1,t2,xmi,ymi)
global ring sid post hsori calmid rm cm gr hsobs xbeg xfin
ybeg yfin obx oby...
xsizob ysizob
hs=hsori;
if(xmi>=obx+1)

if(xmi<=obx+xsizob-1)

88
if(ymi>=oby+1)

if(ymi<=oby+ysizob-1)

hs=hsobs;
end
end
end
end

ring(xmi,ymi)=t0+sqrt((hs‘2)-.25*((t2-t1)‘2));

valuemid=t0+sqrt((hs‘2)-.25*((t2-t1)“2));

B.4 corner.m

function [valuecorn]=corner(t0,t1,t2,xpo,ypo)
global ring sid post hsori calmid rm cm gr hsobs xbeg xfin
ybeg yfin
global obx oby xsizob ysizob inflec cond least
hs=hsori;
if(xpo>=obx+1)
if(xpo<=obx+xsizob-1)
if(ypo>=oby+1)

if(ypo<=oby+ysizob-1)

hs=hsobs;
end
end
end
end

ring(xpo,ypo)=t0+sqrt((2*(hs‘2))-((t2-t1)“2));

89
prevrin=t0+sqrt((2*(h8“2))-((t2-t1)“2));
valuecorn=t0+sqrt((2*(hS‘2))-((t2-t1)“2));
if(gr>1)
for chek=1z10

if(sid==1)

if(least(4*(gr-2)+1,chek)==xpo)
t0=ring(xpo,ypo-1);

t1=ring(xpo+1,ypo-1);
t2=ring(xpo-1,ypo-1);
ring(xpo,ypo)=t0+sqrt(((hs‘2))-.25*((t2-t1)“2));
valuecorn=t0+sqrt(((hs‘2))-.25*((t2-t1)“2));
sidel=1;
xpo=xpo;
ypo=ypo;
Zpause

end

end
if(sid==2)
if(least(4*(gr-2)+2,chek)==ypo)
t0=ring(xpo+1,ypo);

t1=ring(xpo+1,ypo+1);
t2=ring(xpo+1,ypo-1);
ring(xpo,ypo)=t0+sqrt(((hs‘2))-.25*((t2-t1)“2));
valuecorn=t0+sqrt(((hS‘2))-.25*((t2-t1)“2));
side2=2;
xpo=xpo;
YP°=YP°i

Zpause

90

end
end
if(sid==3)
if(least(4*(gr-2)+3,chek)==xpo)
t0=ring(xpo,ypo+1);
t1=ring(xpo+1,ypo+1);
t2=ring(xpo-1,ypo+1);
ring(xpo,ypo)=t0+sqrt(((hs‘2))-.25*((t2-t1)“2));
valuecorn=t0+sqrt<((h8“2))-.25*((t2-t1)"2));
sid93=3;
xpo=xpo;
ypo=ypo;
Zpause
end
end
if(sid==4)
if(least(4*(gr-2)+4,chek)==ypo)
t0=ring(xpo-1,ypo);
t1=ring(xpo-1,ypo+1);
t2=ring(xpo-1,ypo-1);
rinngpo,ypo)=t0+sqrt(((hs‘2))-.25*((t2-t1)“2));
valuecorn=t0+sqrt(((hs“2))-.25*((t2-t1)‘2));
side4=4;
xpo=xpo;
YP°=YP°5
Zpause
end

end

end

end

91

Bibliography

[1] H.D Keith and F..]. Padden. J. Appl. Phys, 35:1270—1286, 1964.

[2] H. D. Keith and Jr. F. J. Padden. Spherulitic crystallization from the melt. ii.
influence of fractionation and impurity segregation on the kinetics of crystalliza-

tion. Journal of Applied Physics, 3521286 — 1296, 1964.
[3] D. C. Bassett and A. S. Vaughan. Polymer, 26:717 — 725, 1985.
[4] M. Avrami. Journal of Chemical Physics, 8, 1940.

[5] Hoffman. J. D. and Lauritzen. J. 1. Journal of Research of the National Bureau

of Standards, 65, 1961.

[6] W. W. Mullins and R. F. Sekerka. Stability of a planar interface during so-
lidification of a dilute binary alloy. Journal of Applied Physics, 352444 — 451,

1964.

[7] W. W. Mullins and R. F. Sekerka. Morphological stability of a particle growing
by diffusion or heat flow. Journal of Applied Physics, 342323 — 329, 1963.

[8] R. Trivedi and W. Kurz. Morphological stability of a planar interface under
rapid solidification conditions. Acta Metallurgica, 34:1663 — 1670, 1986.

[9] SR. Coriell et al. Journal of Crystal Growth, 1402139 — 147, 1994.

92

93
[10] M. E. Glicksman. Fundamentals of Dendritic Growth, volume 1. 1996.

[11] S. Balijepalli and J. M. Schultz. Incipient instability of growth in blends of

polymer melts. Macromolecules, 29:2095—2102, 1996.

[12] A.J.Lovinger and Carl C. Gryte. Model for the Shape of polymer spherulites
formed in a temperature gradient. Journal of Applied Physics, 47:1999 - 2004,

1976.
[13] E. E. Wahlstorm. Optical Crystallography. Wiley, New York, 1969.

[14] Y. Huang and J. Petermann. Spherulite growth of polybutene—l in a thermal

gradient. Polymer Bull, 242649-656, 1990.

[15] T.G.Ryan and P.D.Calvert. Diffusion and annealing in crystallizing polymers.

Polymer, 232877—883, 1982.

[16] Calvert RD. Billingham N .C. and Uzuner A. Studies of diffusion in polymers by
ultra-violet microscopy-i. Eur. Polym. Journal, 25(7):839—845, 1989.

[17] Calvert RD. Billingham NC. and Uzuner A. Studies of diffusion in polymers by

ultra-violet microscopy-ii. Eur. Polym. Journal, 31:258—264, 1990.
[18] D. Turnbull and J. C. Fisher. Journal of Chemical Physics, 17:71, 1949.
[19] R. L. Cormia et al. Journal of Chemical Physics, 37:1333, 1962.

[20] S. Chew et al. The crystallization kinetics of polyethylene under isothermal and

non-isothermal conditions. Polymer, 302874 — 881, 1989.

[21] R. D. Icenogle. Journal of Polymer Science: Polymer Physics Edition, 23:1369
— 1391, 1985.

94
[22] G. Eder, H. Janeschitz-Kriegel, and S. Liedauer. Crystallization processes in
quiescent and moving polymer melts under heat transfer conditions. Prog. Polym.

Sci, 15:629-714, 1990.

[23] F. Bernauer. Cedrillte K ristalle in Forschungen zur Kristalkunde.

[24] Davis G. T. Hoffman, J. D and Lauritzen J. I. Treatise on Solid State Chemistry,

volume 3.

[25] F. P. Price and R. W. Kilb. J. Polym. Sci, 57:395, 1962.

[26] R. Yamada K. Sasaguri and R. S. Stein. Journal of Applied Physics, 3523188,
1964.

[27] Y. Fujiwara. Kolloid Z. Z. Polym., 226:135, 1968

[28] J. M. Crissman and E. Passaglia. J. Res. Natl. Bur. Stand.,Sect. A, 70:225, 1966.

[29] Jin-Young Jung and Michael M. Chen. 1998.

[30] W Kurz and DJ Fischer. Fundamentals of solidification. Trans Tech Publica-

tions, Netherlands, 1992.

[31] RF. Mehl W.A. Johnson. Trans. Am. Inst. Min. Engrs, 135:416, 1939.

[32] AN. Kolmogorov. Izv. Akad. Nauk SSSR Ser. Mat, 1:355, 1937.

[33] Bernhard Wunderlich. Macromolecular Physics, volume 2. Academic Press New

York and London, 1982.

[34] A. Benard. A model for phase change kinetics in polymer composites with sig-

nificant fiber surface nucleation. to appear in Acta Mater., 1998.

95
[35] A. Benard and S. G. Advani. An analytical model for the Spherulitic growth
in fiber reinforced polymers. to appear in Journal of Applied Ploymer Science,

1998.

[36] SD. Hong T. Yuasa R.J. Tabar, A. Wasiak and RS. Stein. J. Polym. Sci. Polym.
Phys. Ed., 19:49, 1981.

[37] E. Piorkowska A. Galeski. J. Polym. Sci. Polym. Phys. Ed, 192731, 1981.

[38] J.M. Haudlin D. Lefebvre N. Billion, C. Magnet. Colloid Polym. Sci, 2722633,
1994.

[39] L. Rebenfeld N.A. Mehl. J. Polym. Sci. 3., 3121677, 1993.

[40] C. Auer G. Hinrichsen Th. Krause, G. Kalinka. J. Appl. Polym. Sci, 51:399,
1994.

[41] Robert C. Mcowen. Partial Differential Equations. Prentice Hall, 1996.

[42] M. Biermann G. E. W. Sculze. Journal of Materials Science Letters 12, 12:11,

1993.

[43] Jaime O. Chua A.J.Lovinger and Carl C. Gryte. Studies on the or and 6 forms
of isotactic polyprOpylene by crystallization in a temperature gradient. Journal

of polymer Science, 152641 -— 656, 1977.

[44] S. Swaminarayan Ch. Charbon. Modeling the microstructural evolution of ther-

moplastic composites. Materials Science and Engineering, 1997.

[45] John Vidale. Bulletin of the Seismological Society of America, 7822062 — 2076,
1988.

[46] J .A.Sethian. Level Set Methods. Cambridge University Press, 1996.

96
[47] B. Engquist and S. J. Osher. Stable and entropy-satisfying approximations for

transonic flow calculations. Math. Camp, 34:45, 1980.

TE UNIV. LIBRRRIES

MICHIGAN srn
[[[llllllllllll[l] [ll
31293 2

CD 037’

Ill]!

11355

II] Ill]
11 ED