Bifurcation and morphological instability


Material Information

Bifurcation and morphological instability
Physical Description:
vii, 120 leaves : ill. ; 28 cm.
Nadarajah, Arunan, 1959-
Publication Date:


Subjects / Keywords:
Bifurcation theory   ( lcsh )
Crystal growth   ( lcsh )
Solidification   ( lcsh )
Crystallization   ( lcsh )
bibliography   ( marcgt )
theses   ( marcgt )
non-fiction   ( marcgt )


Thesis (Ph. D.)--University of Florida, 1988.
Includes bibliographical references.
Statement of Responsibility:
by Arunan Nadarajah.
General Note:
General Note:

Record Information

Source Institution:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
aleph - 001124700
notis - AFM1769
oclc - 20112189
sobekcm - AA00004819_00001
System ID:

Full Text









In the course of the three and a half years of effort that went

into this dissertation, I received help from numerous people in various

ways that contributed to its completion. Enumerating them all would

make this note be at odds with the spirit of conciseness of the rest of

the document and so, reluctantly, I resolved to include only those who

had made a direct contribution. Needless to say, I still remember with

gratitude my debt to the rest.

First mention should be made to my erstwhile mentor, Professor M.S.

Ananth at the Indian Institute of Technology, Madras, who first kindled

an enduring interest in theoretical transport phenomena and encouraged

my proclivities toward graduate study. The greatest help came from my

advisor Dr. R. Narayanan, who apart from getting me involved in morpho-

logical instability and giving many suggestions also imparted to me a

solid background in mathematics and hydrodynamic stability theory, not

to mention helping me in numerous other ways. Professor L.E. Johns,

Jr., first introduced me to linear operator theory and gave suggestions

too regarding my research, but more importantly, he was my "chemical

engineering conscience," broadening my vision when I tended to

specialize too much and keeping the objectives in perspective when I got

wrapped up in abstract theoretical points. It is not an exaggeration to

say that I probably would not have progressed this far academically

without these three individuals.

I would also like to thank Drs. S.R. Coriell and G.B. McFadden of

the National Bureau of Standards for many discussions and suggestions

regarding the subcritical nature of morphological instability; the

members of my supervisory committee: Professor U.H. Kurzweg, Dr. G.K.

Lyberatos, Dr. W.E. Lear, Jr., and Dr. S.A. Svoronos for their time and

effort on my behalf; the Department of Chemical Engineering for provid-

ing a research assistantship during the first two years of my Ph.D. work

and the Department of Mathematics for a lectureship during the last one

and a half. My thanks also go to my uncle Dr. R.S. Perinbanayagam for

being a role model and helping me evolve a "meaningful philosophy of

life" and cope with the stress of graduate school.

Finally, I wish to express my gratitude to my colleague S. Pushpa-

vanam for many "enlightening" discussions and to Debbie Hitt for doing a

superb job of typing the manuscript and "correcting" my Queen's English



ACKNOWLEDGEMENTS .........................................

ABSTRACT .................................................


1 DESCRIPTION OF THE PROBLEM ........................


2.1 Early Work ..................................
2.2 Later Research ..............................
2.3 Inclusion of Other Effects ..................
2.4 Experiments in Morphological Instability ....
2.5 Limitations of Existing Models and
Unaddressed Issues ........................

3 A UNIFORM FORMULATION .............................

3.1 The Formulation .............................
3.2 The Linear Stability Problem ................
3.3 The Adjoint Problem and Exchange of
Stabilities ...............................
3.4 Finite Containers and the Most Dangerous
Wavenumber ................................

4 SUBCRITICAL BIFURCATION ...........................

4.1 Theory .....................................
4.2 The Second Order Problem ....................
4.3 The Third Order Problem .....................
4.4 Calculations and Comparisons ................


5.1 Rayleigh-Marangoni Convection in Brief ......
5.- The Augmented Morphological Problem .........
5.3 Comparison of Morphological Instability
with Rayleigh-Marangoni Convection ........


6.1 Nature of Imperfections .....................
6.2 Imperfection Due to Heat Loss ...............



















6.3 The Outer Expansions ........................ 78
6.4 The Inner Expansions ........................ 83
6.5 Imperfection Due to Advection in the Melt ... 87
6.6 Nonexistence of the Planar State ............ 93
6.7 Asymptotic Solution ......................... 94
6.8 Controlling Imperfections ................... 99

7 NEW DIRECTIONS .................................... 102

7.1 Transition to Dendritic Growth .............. 102
7.2 Extension to Semiconductor Materials ........ 103
7.3 Inclusion of Microscopic Effects ............ 104
7.4 Numerical Methods ........................... 104
7.5 Experiments ................................. 105


A NOMENCLATURE .......................................107

B PHYSICAL PROPERTIES ................................113

REFERENCES ............................................... 115

BIOGRAPHICAL SKETCH ...................................... 120

Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the
Requirements for the Degree of Doctor of Philosophy



Arunan Nadarajah

August, 1988

Chairman: Dr. Ranganathan Narayanan
Major Department: Chemical Engineering

Morphological instability refers to the tendency towards spatial

pattern formation on the liquid-solid interface when a dilute binary

mixture is solidified or fused. The importance of this phenomenon is in

the growth of metal alloy and semiconductor crystals from their melts,

where it influences the solute or dopant concentration resulting in

nonuniform physical and electrical properties.

Previous formulations of morphological instability have involved

several simplifying assumptions which restricted it to the study of a

region immediately surrounding the interface. The models have limited

validity and they require separate treatments for different situations

like freezing and melting. In this study a new uniform approach is

presented which considers the entire melt and crystal domain and is

applicable to all situations. Earlier formulations are shown to be

approximations of this and exchange of stabilities is proven


The model is then used with a weakly nonlinear technique to predict

the shape of the bifurcation diagram for various cell patterns. The

subcritical nature of morphological instability is shown and regions of

its prevalence are determined over the entire domain of experimental

parameters. This was used to compare with experimental results and to

determine optimal crystal growth regions.

A comprehensive comparison of morphological instability with con-

vective instability was undertaken and this phenomenon was shown to

resemble Marangoni convection in its mathematical and physical fea-

tures. This was done in order to introduce some of the multitudinous

mathematical techniques employed in convective instabilities into mor-

phological instability and specifically was used here to complete the

eigenspace of the linearized problem.

Two imperfections which reside in the domain, heat loss at the

container wall and advection in the melt, were considered and shown to

be bifurcation breaking imperfections. Solutions to the problem were

obtained in both cases by matched asymptotic expansions and based on

these results a practical way of minimizing the effect of these imper-

fections was suggested.


Morphological instability refers to the process of spatial pattern

formation at the liquid-solid interface when a binary mixture is

solidified or fused. This is a problem of hydrodynamic instability and

like all other problems of this nature, for this phenomenon too there is

an onset point where the initially planar interface first begins to

deform and forms cellular patterns. These grow into deeper finger-like

shapes and eventually forming side branches and tree-like dendritic


The importance of this phenomenon is in the growth of metal alloy

and semiconductor crystals from their melts. Morphological instability

affects not just crystal shape but also the solute or dopant

concentration, resulting in nonuniform physical and electrical

properties in the crystal. This is especially perfidious in

applications where superfine crystals with very consistent properties

are required. Recently there have been some indications that crystal

quality can be significantly improved by growing them in low gravity as

this reduces other problems associated with crystal growth like natural

convection, but unfortunately not morphological instability. Growing

the crystal at very high temperature gradients or at very low growth

rates avoids morphological instability but crystals grown at high

temperature gradients are of poor quality due to thermal stresses and

the very low growth rates makes the process very expensive. Hence this

becomes a problem not only of avoiding crystal surface deformations but

one of optimization of the process as well.

There are several methods for growing crystals from the melt,

distinguished by the hydrodynamics of growth. The three basic ones are

Bridgman, Czochralski and float zone and most techniques are variations

of these, like horizontal and vertical Bridgman. A typical Bridgman

experiment is shown in Fig. 1-1. The material is usually in a quartz

ampoule and is melted and then recrystallized in the Bridgman furnace.

The upper part of the ampoule is maintained at a higher temperature than

the lower and solidification proceeds upwards and the ampoule is pulled

downwards at the same velocity v. The top of the melt is protected by a

liquid encapsulant like B20 At the end of the process the ampoule is

broken to retrieve the crystal.

In the float zone technique shown in Fig. 1-2, the ampoule or the

material itself is pulled through a circular furnace and the melting and

recrystallization proceeds simultaneously. As this technique can be

done even without a container it avoids the problem of impurities from

the ampoule entering the crystal, but its more difficult to maintain a

uniform temperature gradient. Figure 1-3 shows the Czochralski method

where the crystal is rotated and pulled from a melt pool. As the

emphasis here is on the liquid-solid interface, the modelling of the

crystal growth process will be kept as general as possible but will

resemble the Bridgman technique the most.

The temperature and concentration profiles during typical crystal

growth conditions are shown in Fig. 1-4. The temperature profiles in

the liquid and solid are virtually straight lines and solute

concentration in the solid is virtually constant. But the solute

-Heating Element

-Liquid Encapsulant

-Quartz Ampoule



Pulling Rod

Figure 1-1. Bridgman growth



t Heat

Figure 1-3. Czochralski

Figure 1-2. Float zone

concentration in the liquid CZ changes sharply near the interface

because of solute rejection on solidification. This in turn has an

effect on the freezing point depression in the liquid as shown in the

figure. Figure 1-5 shows the same profiles but in a situation where the

freezing point in the liquid TM now exceeds the actual liquid

temperature (this change can be brought about by either reducing the

liquid temperature gradient or by making the change in C near the

interface even sharper by increasing the growth velocity). This is

referred to as constitutional supercooling and the system responds to

this unstable situation by interface deformation. Countering this is

the interfacial tension which always acts to minimize the surface area

which in this case is the planar surface. When this balance is upset,

or in other words when the onset conditions are exceeded, the interface

loses planarity and forms cellular patterns. Figure 1-6 shows the

profiles for the fusion case, where we could have TM in the solid being

less than the actual solid temperature and once again interfacial

deformation is the system's response, balanced by interfacial tension.

The only difference here is that since solid diffusivities tend to be

several orders of magnitude lower than liquid diffusivities, the solute

concentration profile in the solid near the interface will vary even

more sharply resulting in lower onset conditions for morphological


In the paper of Trivedi and Somboonsuk (1984) there is a series of

photographs from an experiment (see Fig. 1 in their paper) where

succinonitrile/acetone crystals were grown. The first photograph shows

the liquid-solid interface just after onset and has a discernible

cellular pattern. Later ones show the cells becoming deeper, forming

Temperature or


Figure 1-4.

Concentration and temperature profiles during


Temperature or


Ts Solid

Figure 1-5.


C -

Concentration and temperature profiles during

Temperature or

S-TM Solid

Figure 1-6. Concentration and temperature profiles during fusion



C, -

fingers and ending up as dendrites. The arrows in the first photograph

mark the initial perturbations that eventually become dendrites. In

this thesis we will concentrate only on the region near the onset point,

shown by the first two photographs, though dendrites will be mentioned

in discussions.

In the papers by Morris and Winegard (1969) and Tiller and Rutter

(1956) we see another aspect of morphological instability, the variety

of cellular patterns, finger-like shapes, hexagonal cells and variations

of these. Other possible shapes are cylindrical rolls and rectangular

cells but by far the commonly encountered pattern is the hexagonal

one. The choice of the cell pattern is extremely important and factors

determining this choice will be discussed later. Their figures also

show that, unlike other forms of hydrodynamic instability, the number of

cells on a single crystal is in the hundreds.

Though it is customary to model morphological instability in terms

of the temperature and concentration profiles, in reality this

phenomenon seldom exists in isolation; it is usually coupled with fluid

flow in the melt. There are two kinds of flows that occur. The first

is buoyancy driven solutal convection which is caused by the sharp

solute concentration gradients in the melt. The other is the result of

density change during solidification. When solidification occurs there

is a constant rate of volume decrease which causes the melt to move in

to fill the vacated space. This motion is referred to as advection and

the rigid side walls of the ampoule will cause closed streamlined flow

in the melt as a result (see Fig. 6-2). In addition there will be flows

in the melt in Czochralski growth due to rotation and other kinds of

flows in special growth techniques.

Apart from these, several other parameters affect this phenomenon,

the most important of which are due to the fact that most crystals are

faceted; that is, they have a crystal lattice structure. Hence whether

the lattice axis is aligned or not with the growth direction is

extremely important as can be seen from the experiments of Heslot and

Libchaber (1985). Other important considerations are grain boundaries,

wetting of the ampoule wall and the presence of impurities. Also in

rapid solidification, the system will not be at thermodynamic

equilibrium and kinetic undercooling of the melt becomes significant.


2.1. Early Work

The first successful attempt at explaining morphological

instability qualitatively was by Rutter and Chalmers (1953). They

coined the word "constitutional supercooling" to describe the existence

of unstable melt regions near the interface where the freezing

temperature can be higher than the liquid temperature itself and

correctly identified this as the cause of interface deformation.

Tiller, Rutter, Jackson and Chalmers (1953) quantified

constitutional supercooling and for instability came up with the


mG > GC (2.1-1)

where m is the absolute value of the liquidus slope, G the liquid

temperature gradient at the interface and G i the solute concentration

gradient in the liquid at the interface. The negative sign is caused

by Gc and Ge being in opposite directions (see Fig. 1-5). As can be

seen these simple thermodynamic explanations do not take into account

the stabilizing effect of interfacial tension. To do so would require

casting the problem as one of hydrodynamic instability and obtaining the

onset conditions from a linear stability analysis. This is exactly what

Mullins and Sekerka (1963, 1964) did when they considered the problem

with temperature and concentration equations in the liquid and solid and

boundary conditions at the interface. Their criterion for instability

to an infinitesimal disturbance was

mG > Gh + aT X/L (2.1-2)

where GT is the weighted temperature gradient

GT = (k2.GI + kG )/(kz + k ) (2.1-3)

Here ks and k are the solid and liquid thermal conductivities. G- is a

modified liquid concentration gradient given by

G = G c(a v/2D )/(a (1/2 k)v/D ) (2.1-4)

2 2 c 1/2
a2 = (a2 + v /4D ) (2.1-5)

where a is the wavenumber of the disturbance, k the solute distribution

coefficient, TM the melting temperature of the pure solvent, Lh the

latent heat of fusion and X the interfacial tension. This analysis laid

the foundation for all further work in morphological instability.

Following them Woodruff (1968) did the linear stability analysis along

the same lines for the melting problem and came up with the same

criterion as (2.1-2) but with

G = G (na v/2D )/(nka + a. + (1 k)v/2D.)


2 2 2 1/2
a = (a + v /4D2 ) (2.1-7)
5 S

where Gso is the solid concentration gradient at the interface and n is

D /D the ratio of solute diffusivities.

2.2. Later Research

The next major contribution to the problem was made by Wollkind and

Segel (1970) who proved "exchange of stabilities" for this problem for

most parameter ranges. Proving exchange of stabilities is equivalent to

showing the existence of the onset of steady state nonplanar

solutions. They also considered the weakly nonlinear regime after the

onset of instability and using the method of Stuart (1960) and Watson

(1960) analyzed the problem for the case of two-dimensional rolls,

showing the existence of subcritical bifurcation at most growth


The method of Stuart and Watson is essentially the theory of Landau

(see Drazin and Reid (1981)) and follows the dominant mode of

instability into the weakly nonlinear regime. In this form it is

applicable only to disturbances of one cellular pattern at a time, but

Segel and Stuart (1962) extended this theory to the prediction of the

preferred pattern by considering the interaction of two specified modes

of disturbance. Depending on the way these two modes were combined it

was possible to obtain two dimensional rolls or hexagonal cells and they

showed that the experimental parameters would determine the stability of

these patterns. Sriranganathan, Wollkind and Oulton (1983) adopted this

method for morphological instability and gave parameter ranges where

each type of cell was stable. The limitation of the method is that it

considers hexagonal and two dimensional roll patterns but not

rectangular cells or cylindrical rolls and ignores the effect of

container shape and size which have been shown to be important in wave

pattern selection (see Koschmeider (1967)).

Ungar and Brown (1984a) considered the highly nonlinear problem and

after making several simplifications obtained solutions using the finite

element method. Finite elements can handle highly nonlinear problems

and give very accurate numerical solutions but are extremely time

consuming. Solving the full morphological problem is a very expensive

proposition by this method and hence Ungar and Brown simplified the

problem by ignoring the latent heat and solid diffusivity and assuming

that thermal conductivities in liquid and solid were equal. This

allowed them to reduce the problem to a "one-sided model" consisting of

variables in the liquid region only, considerably simplifying the

algebra and saving computer time. Such a model will have only limited

validity in highly nonlinear regions and this was borne out when Ungar,

Bennett and Brown (1985) solved the complete problem. But their most

extensive calculations were done only for the one sided model and hence

this is of chief interest. These were done only for the case of two

dimensional roll disturbances and here they showed that contrary to that

reported by Wollkind and Segel, there were multiple regions of

supercritical and subcritical bifurcation. More importantly they showed

that at large deformations of the interface secondary bifurcations

occurred. Ungar and Brown (1985) also modelled the formation of deep

cells in an attempt to follow the transition to dendritic growth.

Nonlinear finite difference calculations were done in a more

limited way by McFadden and Coriell (1984) for the two dimensional

case. Later McFadden, Boisvert and Sekerka (1987) extended the

calculations for the three dimensional patterns of hexagons and cross

rolls. In both cases the enormous expenses involved restricted

calculations to a few parameter values.

2.3. Inclusion of Other Effects

While these workers were investigating the basic problem others

were busy trying to incorporate various influences. The most important

concern was the effect of fluid flow. Delves (1968, 1971 and 1974) in

attempting to approximate the influence of advection and stirring in the

melt, studied the influence of plane Couette flow on the problem. He

showed that two dimensional roll disturbances in the flow direction were

stabilized but there was no effect on disturbances perpendicular to the

flow. Coriell, McFadden, Boisvert and Sekerka (1984) modelled Couette

flow more systematically and came to the same conclusion. Recently

McFadden, Coriell and Alexander (1988) examined the effect of plane

stagnation flow on two dimensional disturbances perpendicular to the

flow and here too the flow as found to be stabilizing.

In another very important development Coriell, Cordes, Boettinger

and Sekerka (1980) studied morphological instability with solutal

convection. They showed that the two in-.;abilities were essentially

decoupled with the melt being unstable to convective disturbances of

long wavelengths and the interface unstable to nonplanar disturbances of

small wavelengths. Also at low growth rates the dominant instability

was convective and the interface was not easily disturbed. At high

growth rates the roles were reversed and at an intermediate velocity the

two instabilities became comparable. It was only at this rate the two

instabilities interacted and the result was the prevalence of

oscillatory instabilities. Their conclusion was that, except at this

particular growth rate, it is usually sufficient to study only the

dominant instability near its onset.

Following Coriell et al. several workers have looked at special

aspects of these two instabilities and their work has been reviewed by

Glicksman, Coriell and McFadden (1986). They all confirmed or refined

the work of Coriell et al. but all the main conclusions mentioned above

still hold.

Several other influences apart from fluid flow have been

incorporated into the model but only a few relevant ones will be

considered here. Coriell and Sekerka (1972, 1973) tried to include the

effect of grain boundaries on morphological instability by assuming that

its only effect was to shift the onset conditions. They failed to

observe that in the presence of grain boundaries there could be no

planar solutions to the problem and that the interface will be nonplanar

at all times. Ungar and Brown (1984b) obtained the solutions to this

problem by matched asymptotic expansions for small grain angles and

using finite elements for solutions of large grain angles.

In rapid solidification kinetic undercooling of the melt is

significant and Seidensticker (1967) included this and showed that it

caused a shift in the onset conditions. The significance of this was

shown by Hardy and Coriell (1968, 1969 and 1970) when they observed

morphological instability in the growth of ice crystals. Constitutional

supercooling was not a factor here and it was shown that kinetic

undercooling was the primary cause. This dual cause for morphological

instability is somewhat analogous to the situation in natural convection

where we find that the variations of density and surface tension with

temperature can both cause convective instability.

2.4. Experiments in Morphological Instability

The early work on modelling morphological instability was prompted

by experimental observations but beyond that very few quantitative

experiments have been done near the onset conditions. This is an

unfortunate state of affairs and experimental verifications of

theoretical predictions are badly needed if further concrete progress on

the theoretical front is to be made. The work of Morris and Winegard

(1969), Trivedi and Somboonsuk (1984) and of Heslot and Libchaber (1985)

have already been mentioned. Recently de Cheveigne, Guthman and Lebrun

(1985, 1986) have attempted to verify the weakly nonlinear and strongly

nonlinear theoretical predictions and one hopes that more work along

these lines will follow.

De Cheveigne et al. performed their experiments on

succinonitrile/acetone and CBr4/Br2 systems. (These organic mixtures

are much easier to work with than metal alloys as they are generally

nonfaceted, transparent and require small temperature gradients and

hence they have been very popular with experimentalists.) They found

that the cell pattern formed and its dimensions were strongly dependent

on geometry of the container. More importantly when they ran the

experiments for two dimensional roll patterns, they found only

subcritical instability.

2.5. Limitations of Existing Models and Unaddressed Issues

In Chapter 1 the cause of constitutional supercooling was explained

as being due to the sharp solute concentration gradient in the liquid

near the interface, while elsewhere in the liquid and the solid the

solute concentration was practically a constant. It would seem then

that the only region of interest is the interface and a liquid "boundary

layer" adjacent to it. This has prompted all previous workers in

morphological instability to consider D /v as the characteristic length

of the problem and to ignore solid diffusion. A typical value

of D /v is 100 microns and this means that the far ends of the melt and

crystal are infinitely far away and the domain of the problem is

effectively confined to the liquid boundary layer mentioned above. For

the melting problem a characteristic length of Ds/v is used and the

domain becomes an even smaller boundary layer in the solid.

These assumptions considerably simplify the algebra involved and

hence their popularity. But they constrain the validity of the model in

several ways. The most obvious one is that they necessitate the melting

and solidification problems to be studied separately, even though they

only differ in the direction of the growth velocity. Besides this

assumption fails for very small growth velocities, as it introduces a

singularity at v = 0. Later we will show that neglecting solid

diffusion also introduces a singularity and makes the model fail in the

nonlinear regime.

Finally, any effect which resides in the entire domain, not merely

the boundary layer, cannot easily be incorporated into the model, which

is why all influences on morphological instability studied so far are

either boundary layer effects (e.g., solutal convection) or interfacial

effects (e.g., grain boundaries and kinetic undercooling). Phenomena

that span the entire domain, like advection in the melt or imperfect

insulation of the ampoule walls, have been either inadequately treated or

ignored completely. Hence there is a need for a model that includes the

entire liquid and solid domains which would be applicable for all growth

velocities. This model should also dispense with the separate

treatments accorded so far to the solidification and fusion problems

with one uniform formulation.

In Section 2.2 it was mentioned in connection with the work of

Wollkind and Segel (1970) and of Ungar and Brown (1984a), that this

problem oscillates between subcritical and supercritical instabilities

for the case of two dimensional roll disturbances. They did not,

however, compute the ranges of each type of instability for the

experimental parameters involved. This is necessary in light of the

experiments of de Cheveigne et al. (1986) who observed no supercritical

instability. Also the extension of these predictions to three

dimensional disturbances like hexagonal and rectangular cells is yet to

be done.

It would not be an exaggeration to state that the inspiration for

all the theoretical work done so far in morphological instability has

come from Rayleigh-Benard convective instability. A comparison between

the two problems would be invaluable as a source of continued

inspiration and as a way to draw conclusions and conjectures about

morphological instability from the vast published literature on

Rayleigh-Benard convection. Hurle (1985) has attempted this but his

work can only be regarded as perfunctory and there exists a need for a

more rigorous treatment of the issue.


3.1. The Formulation

Since we are not making the assumption that the liquid and solid

are very deep, the problem has to be formulated very carefully,

especially with regard to the outer boundaries, if we are to avoid an

intractable moving boundary problem.

A typical crystal growth set up is shown in Fig. 3-1. The ampoule

is heated by the heating coils surrounding it and they keep the melt

region at a temperature T1 and the crystal at T2. The temperatures T1

and T2 are maintained constant by means of thermocouples located at z =

s and z = -X. The ampoule is pulled towards the cooler end at the same

velocity V at which the crystal grows, thus keeping the interface

stationary. The region near the interface is protected by an insulating

shield and it is this region that becomes the domain in our model.

So in this model the outer liquid and solid boundaries become fixed

at z=-1 and z=s respectively and the solid will be moving with a bulk

velocity v and the liquid with a bulk velocity v/Y, where Y is the ratio

of densities p /P In this section we will assume that Y=1 and

consider the effects of Y not being unity in Chapter 6 as this would

cause advection and fundamentally alter the basic problem. Also we will

assume that the melt concentration at the outer liquid boundary is a

constant C1.


z = 0

So o o o o 0 0 0 0 0 0

1~~ v //1 77 7 7777 -7 7f/. / f/ If/./If///7777


z = -I


o o o o o o o o o o o o

z= S

Heating Coil

Figure 3-1. Experimental set up



The domain equations in the liquid melt are



az x

3z 2 z
T at. De \



In the solid region the equations are



aT 2

ac 2
v = D C
5z s s

where T and C are the temperature and solute concentration, D

diffusivity and thermal diffusivity, with the subscripts

referring to the liquid and solid.

The boundary conditions are

Tz = T1, C C at z = -1

T =T at z=s
s 2



and a the

II and s



At the liquid-solid interface we will use r to denote the departure

from planarity and write the boundary conditions at z-=

T = T = T mC -T -H (3.1-7)
s M I M Lh

k VT n kT n = L(v + D ) (3.1-8)
2. 9. Ss h +Dt)(.18

kC = Cs (3.1-9)

D VC* n (v + )C D VC n (v + )C (3.1-10)

where TM is the melting point of the pure solvent, A the interfacial

tension, Lh the latent heat, m the absolute value of the liquidus

slope, ki and ks the liquid and solid thermal conductivities, k the

distribution coefficient, n the normal at the interface directed into

the solid and H is the curvature of the liquid-solid interface and in

Cartesian co-ordinates is given by

2 2 2
H [ 'sC (1 + (i)2) 2 S + 1 (1 + (LC)2 )]
S2 ax ax ay axay 2 ax
dx yy

[1 + 2 + (C))2 2-3/2 (3.1-11)
ax Dy

In cylindrical co-ordinates, if we assume circular symmetry, it


H- [ 12 (1 + (,)2) a[1 + (,)2]-3/2 (3.1-12)
2 Dar rar 3r

It is assumed that the side walls are sufficiently far apart and

well insulated to enable us to impose periodic boundary conditions in

that direction.

To convert these equations into the dimensionless form we will use

the liquid depth Z as the characteristic length and the diffusive time

D /v2 as the characteristic time.

The temperature will be made dimensionless by T = (T-TM)/G t with

GT = (kG + ksG )/(ks+ k ) (3.1-13)

where G, and Gs are the temperature gradients in the liquid and solid

in the quiescent planar state.

Similiarly the concentration will be made dimensionless by

C = C/G a with

G = (D G + DG )/(D + D ) (31-14)

where G.o and Gso are the concentration gradients in the liquid and

solid in the quiescent planar state.

So in dimensionless form, if we neglect the Lewis numbers D /a in

the temperature equations, we have

2i =0 (3.1-15)

ac aC (3.1-16)
+ v = C
at 3z

V2T =0 (3.1-17)

DC aC (3.1-18)
+v r ni EC
5t zs s

where v and = D /D

The boundary conditions are

Ti = (T TM)/GTL = T1


C = 1 /GZ = C
i 1 C

at z=3s/=s,

Ts = (T2 T )/G I T2

at the interface z=

Ti = = -SeC -A f

VT n BVT n = L( + )

kCi = s

VC n (v + )C = VC n (v + )
Dt Dt






where is the total derivative, a the ratio of thermal conductivities
k /ki and Se, A and L are the Sekerka, capillary and latent heat numbers


Se = (3.1-25)

A M 2 (3.1-26)

L =k I (3.1-27)
kiGT z

at z= -1,

At this point a discussion of the choice of a critical parameter

becomes imperative. The experimentally variable parameters for this

problem are Gt, C1 and v or their equivalents in this formulation

GT, Go and v. Most previous workers have adopted G or C as their

critical parameter but Hurle (1985) has proposed Se, by analogy with the

Rayleigh-Marangoni problem. Recently de Cheveigne et al. (1986) have

advocated the use of v from an experimentalist's perspective. Although

this is a valid choice the reason no one else has used it so far is

probably because v occurs in the domain equations and will give rise to

an infinite number of eigenvalues in the linearized problem. We second

Hurle's suggestion and choose Se as this seems to be the naturally

occurring coupling factor between temperature and concentration and the

fundamental cause for morphological instability. Besides it includes

both GT and Go But the suggestion of de Cheveigne et al. still

remains a valid one.

3.2. The Linear Stability Problem

From now onwards the for dimensionless variables will be

dropped. The steady state planar solution to this problem occurs when

c (x,y) or c (r) = 0 (3.2-1)

The solution is then

T c(z) = (T1 + SeC2 (o))z SeC2 (o)


T (z) = (T2 + SeCc (o))z/s SeC c(o)

C (z)= C [ (1-k)exp(vz) 1]/ (1-k)exp(-v) + 1]
to 1k-exp(-v(1+s/n)) k-exp(-v(1+s/n))

C (z) C r (1-k)exp(vz/n) + (1-k)exp(-v) +
so1 kexp(v(1+s/n))- 1 k-exp(-v(1+s/n))




To write the equations of the linear stability problem we will

impose an infinitesimal disturbance on the steady state solution.

T = Tto + T
2. ic 9&


with similar expansions for Ts, CV, Cs and .

Considering the linear stability problem, in order to separate

variables we will assume a horizontal cell pattern. This pattern for

two dimensional rolls is given by

2irnx 2wn
c (x) = Cos with wave number a =
n L n L

2inx 2irny 2irn
for cross rolls 0 (x,y) = Cos -L- + Cos a
n L L n L



for rectangular cells

0n(x,y) = Cos 2Ln Cos 2 an = 2rn(1/L2 + 1/L2) 1
n L1 L2 n 2


for hexagonal cells

2wnx 2iny + ny C nay
r (x,y) = 2Cos 23L Cos 3L + Cos L, an 4
n 73-L 3L 3L -n -3L

for cylindrical rolls n (r) = J (a nr/R)
n o n



where an are the zeros of J1, and Jo and J1 are the zeroth order and

first order Bessel's functions of the first kind. We can now write

T (xy,z,t) = T1V(z) 41 (x,y)eot


with corresponding forms for the other variables, where a is the

eigenvalue of the linear stability problem.

The linear stability problem becomes in the domain

2 2
(D a )T = 0

2 2
oC2= (D vD a )C

2 2
(D a )T = 0

- C





2 vD 2 -
=(D --a )Csl

Here we have used D to denote -
The boundary conditions are

T =0, C1 = 0

T = 0

at z = -1


at z = s (3.2-18)

At the interface, z=0

T Ts + ;1 (G s) = 0

T + SeC + G (C + SeG c- aA ) =

DTl -O DTsl = aL 1

K1 C + C 1 (kG c Gso) = 0

DCiM nDCs1 vC 1 + vCl1 = Gv(C 0- Cso c01






Solving the system we obtain a general equation for morphological


Ce + (1-k)tanha stanha
Ce k(na + v/2 tanha s)tanha + (a -v/2 tanha )tanha s
S S Xf C X* s

= G + aA -

a(k tanhas + k tanha)
A s

a a a2 [2 v2 1/2
s n 2

a, = [a2 + + v2/4] 1/2

k G tanhas + k G tanha


k tanhas + k tanha
II s






G c(a- v/2 tanha )tanha S + G (na v/2 tanha s)tanhaZ
c (a v/2tanha )tanha s + k(na + v/2tanha s)tanhaZ
R A 3 a S &

If we can show exchange of stabilities for this problem, then for

neutral stability.

GT + a2 A
Se = (3.2-29)

Here we have already used the fact that the Se obtained from the neutral

stability curve will be the same as Seo defined later in eqn. (4.2-2).

In eqn. (3.2-29) if we let v become very large the critical wave number

amin (for which Se0 is a minimum) also becomes very large and we can
approximate all the tanh terms to unity. Further if we also neglect

solid diffusion, n and G are zero and the equation reduces to the well

known results obtained by Mullins and Sekerka (1964).

G (a v) oC. (1-k) 2
Seo 2 ] = 1 + a (3.2-30)
a + v(k ) a + v(k ) a(k + k)

For the case of neutral stability this becomes

(1 + a2A )(a /v + k I)
Se = (3.2-31)
G(a/V -2 )

It must be pointed out that setting all the tanh terms to unity is

equivalent to using the diffusive length as the characteristic length.

This is a boundary layer approximation not unlike that used in the study

of pipe flow at high Reynold's numbers.

In our derivations we did not restrict v to be positive and hence

eqn. (3.2-24) is also valid for negative velocities, that is the fusion

problem. If we replace v with -v and here too assume that v is very

large and set the tanh terms to unity we obtain the result of Woodruff


G c (a -1/2) C (1-k) + -
nka + a + v(1-k) nka + a + v(1-k) a(k + k )


and for a = 0

(1 + a A)(nka + a + v(1-k))
Se = ) 2 (3.2-33)
G (na v)
so a 2

To compare this model with the approximations of Mullins and

Sekerka and that of Woodruff, Seomin was calculated for various growth

velocities for the Pb-Sn system and the results are shown in Fig. 3-2.

As can be seen the approximations hold up very well for most growth

velocities but begin to fail for small velocities. (The thermophysical

data for the Pb-Sn system were those of Coriell et al. (1980).) The

ratio n for the Pb-Sn system is of the order 10-5 but for systems with

much smaller values of n like the C-Austenite system (see Clyne and Kurz

(1981) and Wolf, Clyne and Kurz (1982)) the approximations begin to fail

at higher velocities.

ig U) I s

0 >
4 8


-4 0

o o C

qo 0C
0 0

a -o

o o
?* \ Io -

a 0
O > a

qt r-'-
W ki

3.3. The Adjoint Problem and Exchange of Stabilities

Exchange of stabilities refers to the nonexistence of time periodic

infinitesimal perturbations. Time-dependent infinitesimal perturbations

about the planar state will generally have periodic and nonperiodic

components, with a = or + iai. For some problems it is possible to show

that a. is zero and this is called exchange of stabilities (see looss

and Joseph (1980) for details).

We still have to prove exchange of stabilities for this problem but

before we can do that it is necessary to obtain the adjoint problem. To

accomplish this we will define a column vector 9 and a matrix operator




Se 2 2
-- (D vD-a a)


2 2
(D -a )2


Se 2 2
--(nD -na -oa-vD)


and an inner product < *> = (T T + CC )dz + (TT + CC )dz
0 io ss se
(1 3

L =

2 2
(D a )


where T denotes the complex conjugate of T, T* the adjoint function of

T and

x = (kNG0- Gsc)/(Ge- G5)


Then the domain equations can be written as L*1 = 0 (3.3-5)

In this inner product the adjoint problem becomes in the domain

2 2 ^"*
(D -a )T

2 2 ^*
(D -a )T = 0

2 2 ^
(D + vD a )Czi

2 vD 2 ^A*
in sl

a *
= aC

= c
n s1





Subject to the boundary conditions

T1 C 1
R1 L1

T = 0

at z= -1

at z=s



At the interface z=0

BT1 = T


nCt1 = c1 (3.3-13)

kDCs1 DCl + (kG to- Gs) = 0 (3.3-14)

DT+l D1 + (G G) = 0 (3.3-15)

DT + SeDC + (Gt + SeG aA) = T
51 Si (G SeG (G,- G5) 11

GCCt (1-k)
-(kG SeC (3.3-16)

So far we have been unable to show exchange of stabilities for this

problem directly. But Wollkind and Segel (1970) have proved it when the

boundary layer approximation is valid and we will prove exchange of

stabilities by performing an asymptotic analysis around the boundary

layer solution, which corresponds to n=0 and Pe=0, where Pe is the

Peclet number given by D ./v.

^ 00 ^10 01
T1 = T + nT1 + PeT1 + .... (3.3-17)

00 10 01
a = o + no + Pea + .... (3.3-18)

The other variables are expanded similarly.

If we use L to designate the operator L in (3.3-2) when n and Pe are

zero, then the linear perturbed systems become

L 10 =10 & L 00 f01 (3.3-19) & (3.3-20)
1 1


With boundary conditions

T10 ^10


(3.3-21) & (3.3-22)

at z=0 the equations are

00 ^10 ;10 10
B (0 ) = h

00 ^10 ^01
& B (0 1 1 l )

= 0 1

(3.3-25) & (3.3-26)

where B00 is the boundary operator defined

by eqns. (3.2-19) (3.2.23)

when n and Pe are zero and

^00 10
- (G -
^00 10
- (G +
10 "00
a LS

SeG )

^00 10 10
- N(kG G )
10 00 ^00
a C (1-k)1 +
to 1

00 10 ^ 00
a C (1-k)
toc 1



& f01

01 00

at z = -


at z=-


h =


h will have a similar form with 01 superscripts replacing the 10. It

00* ^00*
we let L and 0 represent the adjoint operator and adjoint
00 ^00
function of L and #* respectively, we can use the solvability

condition on the above system.

00 00 10 ^00* 10
< L > = < f > (3.3-28)
1 1 1

00*^ 00* ^10
=0 (3.3-29)

Subtracting we get

^00* ^00* 10 <00* 01
J(41, ^\ h 0) = <0 f > (3.3-30)

where the lefthand side of eqn (3.3-30) is the bilinear concomitant

evaluated at z=0.

^00* ^00* 01)=0 <00 01
Similiarly J(01 h 0) = <* f > (3.3-31)

^00* ^00
We note that *, and ^o are real as they correspond to a state where
10 01 10 01
exchange of stabilities has been proved, while h10 h f and f are
10 01
also made up of real quantities except for o and o Hence we

conclude from eqn (3.3-30) that a is real and from (3.3-31)

that a is real; i.e. exchange of stabilities holds for the generalized

morphological instability problem upto 0(n) and 0(Pe). So we are

justified in using the neutral stability curve (3.2-29) to calculate Seo

at least up to such order even though we often extend it further. Even

when the boundary layer approximation is valid, exchange of stabilities

does not hold for all growth conditions (see Coriell and Sekerka (1983))

and care must be exercised when such extensions are made.

3.4. Finite Containers and the Most Dangerous Wavenumber

Appealing to the proof of Section 3.3, henceforth we shall only

consider steady state solutions. A typical Seo versus a diagram is

shown in Fig. 3-3 and Seomin is the minimum value of Seo and the

wavenumber at which this occurs is amin. If we can maintain the growth

conditions such that Se < Seomin we can at least say that the planar

interface is stable to infinitesimal perturbations.

As Seomin is the least value of Se at which the planar solution

loses the stability, amin is the wavenumber of the disturbance which is

most likely to occur. Hence this wavenumber is commonly regarded as the

most dangerous and is the wavenumber at which morphological instability

is usually studied. In the derivations and discussions that follow this

is the wavenumber employed and it becomes our "operating wavenumber."

The wavenumber is a 2w multiple of the reciprocal of the

wavelength. If the domain being considered is regarded as being finite

with periodic lateral boundary conditions, the wavelength R (or L

depending on the wave pattern used) in its dimensionless form is now the

aspect ratio, the ratio of the ampoule radius to the melt depth. For

this situation Seo takes on different values depending on the number of

cells formed on the crystal surface as shown in Fig. 3-4 (see also

Rosenblat, Homsy and Davis (1982)). For each value of R there is a

fixed number of cells which is the pattern that is most easily

disturbed, except at certain values of R where two different patterns



Seo min


Figure 3-3. Se vs. wavenumber diagram



Seo min

Figure 3-4. Se curves for various aspect ratios




are equally dangerous. These "horizontal multiple points" are very

important and will be discussed later. At other values of R the most

dangerous number of cells is denoted by N and the corresponding

wavenumber of each cell and the value of Seo are denoted as aN and SeoN

respectively. For the analysis in Chapter 6 where the effect of a

finite container on morphological instability is considered, R is chosen

such that SeoN is as close as possible to Seomin so that the worst

possible case can be examined.


4.1. Theory

The linear stability analysis will only give the onset conditions

for morphological instability. Nonlinear calculations are necessary to

determine the behavior beyond this point. Ideally numerical

calculations should give the most amount of information by being

applicable for small and large deformations, but as can be seen from the

work of Ungar and Brown these calculations are very expensive to carry

out and they were forced to simplify the nonlinear problem and perform

calculations for very few experimental conditions. McFadden et al.

(1987), even though they did not attempt calculations for larger

deformations, were faced with the same restrictions.

This then is the case for weakly-nonlinear methods. They are

generally valid only in a small region very close to the onset

conditions but they can be used to predict the shape of the nonlinear

curve for larger deformations and for several applications this

information is sufficient. More importantly due to the analytical

nature of the techniques they can be used to predict the weakly

nonlinear behavior for all experimental conditions. A case in point is

the work of McFadden et al., most of whose predictions could have been

obtained more cheaply for all parameter values from weakly-nonlinear


Probably the most useful information generated by these theories is

the subcritical behavior of the nonlinear curve. Some typical

bifurcation diagrams are shown in Figs. 4-1 to 4-4. The e = 0 axis

corresponds to the planar solution and initially for small values of Se

the planar solution is stable and usually the only possible solution.

For Se Z SeoN the planar solution becomes unstable and a nonplanar

solution bifurcatess" from the planar one. Figure 4-1 shows a symmetric

bifurcation diagram where nonplanar solutions do not exist for Se <

SeoN. For Se > SeoN, even an infinitesimal perturbation will make the

solution jump from the unstable planar solution to the stable planar one

while for Se < SeoN the planar solution is stable to all

perturbations. This behavior is referred to as supercritical

bifurcation and for these curves Se, = 0, Se2 > 0 (where Se, = dSe/de
2 2
and Se2 = d Se/de at e= 0). Figure 4-2 is nonsymmetric and as can be

seen stable and unstable nonplanar solutions exist for Se < SeoN, making

the planar solution stable to infinitesimal perturbations in this

region, while a large perturbation can make it jump to the stable

nonplanar branch. This is called a subcritical bifurcation diagram and

is characterized by Sel 0. In this situation it is obvious that

growing the crystal at Se < SeoN is no guarantee of avoiding

morphological instability.

Figures 4-1 and 4-2 display the behavior usually seen in most

problems of hydrodynamic instability. Morphological instability is

unusual in having nonlinear curves shown by Figs. 4-3 and 4-4 as well.

These curves have been labelled "backward bending" to distinguish them

from the usual "forward bending" curves. (Actually they are Janus-like

in appearance bending backwards and forwards.) From the point of view

- Stable
--- Unstable


Sei = 0
Se2 > 0

Figure 4-1.

Forward bending, locally
symmetric, supercritical
bifurcation diagram

Figure 4-2.

Forward bending,
unsymmetric sub-
critical bifur-
cation diagram

S- Stable
--- Unstable

---; ----S-



Sei = 0
Se2 > 0

Figure 4-3.

Backward bending, locally
symmetric, subcritical
bifurcation diagram

Figure 4-4.

Se1 # 0

Backward bending,
unsymmetric, sub-
critical bifurca-
tion diagram


Sel # 0

of the crystal grower this is unfortunate as their subcritical nature

increases the occurrence of subcritical bifurcation. Figure 4-3 is the

symmetric case with Se1 0 and Se2 < 0 and Fig. 4-4 the nonsymmetric

one with Se 0.

The symmetry or nonsymmetry of the bifurcation diagram in

hydrodynamic stability is dependent on the cell pattern (see Joseph

(1976)). For morphological instability it will be shown that two

dimensional rolls, rectangular cells and cross rolls produce locally-

symmetric bifurcation diagrams while cylindrical rolls and hexagonal

cells produce nonsymmetric diagrams. It will also be shown that both

backward bending and forward bending will occur for all cellular

patterns depending on the experimental parameters used. Hence

bifurcation can be subcritical or supercritical for two dimensional

rolls, rectangular cells and cross rolls but for hexagonal cells and

cylindrical cells bifurcation is always subcritical.

4.2. The Second Order Problem

Here we begin our weakly nonlinear analysis and in this section we

will calculate Sel, the first derivative of Se with respect to E.

Considering the neutrally stable nonplanar solution near the bifurcation

point Seo, we will expand the variables around the planar steady state


T T i + ET o+ E2T 3 ...... ........... (4.2-1)

Se = Se + eSe + E2Se2 + ................ (4.2-2)
0 <

where e = < 1- /> 2
c c


We have already obtained the linear perturbed solution in Section

3.2. Substituting these expansions into the steady state versions of

eqns. (3.1-15) (3.1-24) and collecting the terms of order e2 we get

the second order perturbed problem.

If the first order perturbed variables are written, following eqn.

(3.2-12), as

T o(x,y,z,Se ) = T 01(zSe o (x,y))


Then the solution to the second order problem becomes

Tv = E Ttn(zSel, T O) )n(xy) (4.2-5)

Substituting (4.2-4) and (4.2-5) into the second order problem and

taking Fourier transforms horizontally, we get

L11 0


and T11 Call = 0

Till =0

B(1 C ) =h11

at z = -1

at z =s




at z=0

where L and B are the same as that used in (3.3-2) and (3.3-25) but with

=0O and Se = Se and

o01(DT01 DTSo1 11
Coi(DTo + Se DC ) vG Se 2 o ]I Se (CO + C0 G
01 01 o 101' 2 to 0 01 11 1 +01 01 to

-a 01 01- T01 11 + 01(T101- BT01)12

-C01( DCo DC ) v12 (kG G / n)
01 101 Sol 2 01 to so 11
-a2c01 (C1l nC )1 + 01 (C01 -C s)I12

L 2
2 f
o o
11 L= L

0 0
o o

3 (x,y)dxdy

2 (x,y)dxdy

"2 "1 S 1 )2 + ( 2]dxdy
fx ay
12 o o (4.2-12)
12 L2 L
2 f (x,y)dxdy

It can easily be seen that for two dimensional rolls, cross rolls

and rectangular cells that Ill = 112 = 0. Hence from the solvability

condition for these patterns Se1 will be zero; that is, the bifurcation

will be locally symmetric. But for hexagonal cells and cylindrical

rolls I,, and I12 will be non-zero and hence Se1 will also be nonzero

and we have nonsymmetric bifurcation and the existence of subcritical




We will analyze the case of Se 0 by considering hexagonal cells

as an example, but the results obtained are applicable to cylindrical

rolls as well. It should be noted here that all the horizontal cell

patterns we have been using are possible solutions to Helmholz's


12 a2 + a 2 = 0 (4.2-13)
2 2
ax ay

and the solutions always come in pairs. For hexagonal cells the

complementary solution to 11 given in (3.2-10) is i

1 = 2 Cos 3-x Sin -Jy Sin 41y (4.2-14)

So the general form of eqn. (5.4) will be

Tto T01( + Pi1) (4.2-15)

To determine p the procedure outlined by Joseph (1976) will be

used. We will proceed the same way as above but using (4.2-15) instead

of (4.2-4) and multiplying by (1 and integrating horizontally. This

will give the same set of eqns. (4.2-6) (4.2-10) but with

L /3L 2
jf J (0 + p1) 21 dxdy
S= 0 (1-p)2 (4.2-16)

f ( + p1 1dxdy
0 0

,/3L 3x 3x

S(0 1 + pi1)o1dxdy

(1 + p "1 2]0 dxdy
ay By 1 a 2
---------= 2- (1p


which can be used to calculate Se1.

We then repeat the process

with *1 and equate the two Se l's obtained. This will result in a cubic

equation for p.

p3 3p 0 (4.2-18)

Se = Se (1-p2)


Se1 G (k=G. G ) (a + v/2tanha )B (na v/2tanha s)
01 (k + ) tanha z ntanha s 11

(1+nB)tanha s
(na + v/2tanha s) 12' 1


B(G G ) a
h s r 11
1 (B+A) tanha

a(GC + Seo G aA)(1/-1)
(G -G )tanha I11
X> 3

(1+A)tanhas I + h
(a+A) 12 2

Se I a(8-1)(kG G)
h [ -o 11 [o--s (BG + G /n)]
2 (k+B) (B+A)tanha 2c sc

12 IL




where B = (aI v/2tanha )tanha s/(na + v/2tanha s)tanhaZ (4.2-23)

and A tanhas/tanha (4.2-24)

There are three solutions 0, /3 and -/3 for p, but it can be seen

that Se1p for /3 and -/3 coincide. It is also obvious that while any

two of the solutions are independent the third is a linear combination

of the other two. Here is an instance where the problem exhibits a

multiplicity of solutions for the same eigenvalue, and could cause

secondary bifurcations further along the bifurcation curve.

The next obvious question is to ask if there is any point at which

Sel in (4.2-20) goes to zero and hence causing Selp to go to zero.

Calculations done for the C-Austenite (i.e. steel) system did give such

a curve (see Fig. 4-5) though for most growth conditions this curve lies

far away from the critical (or most dangerous) wavenumber amin'

intersecting it only at high velocities.

4.3. The Third Order Problem

In Section 4.2, it was shown that for 2-D rolls, cross rolls and

rectangular cells Sel was zero, which implied a symmetric bifurcation

diagram. But to learn more about the nature of this diagram it is

necessary to go to the next order. If we repeat the procedure described

for the second order problem for the terms of order E3 we will obtain

the third order problem.


L 21 = 0

20 40
VELOCITY (pms-1)

Figure 4-5.

Growth velocity vs. wavenumber chart for
C-Austenite system for hexagonal cells






T121 9 C21 0


B(21 21) = h 21

at z= -1

at z=s



at z=0

a 2T -
2 01 2 01

21 2 01 Z01

01 )121

+ SeoC + 2
0 9.01 201Se0

- (Gtot 101Gei)Se2 -

DC 1- 2 3 SeoGC1
o01/ 6 01 0oI c
A3 (3/2 a2 1- I )
01 22 23

- Cs01 + (k DC01 DC01n)

+ v1;(kGc G /n2)
6 01 io se

- 22

- DCs01 21

O1 (DCo01

- nDC01 22
s01l 22



f i dxdy
o o -
21 L2 L1

fI f 22dxdy
o o

L2 L1

0 0

1 ax


+ (- 2dxdy


0 0 dxdy
o o



-I[21 a 02(kC01

'( 2 ) o2 -
By ax y2
ay ax a7

2 21

a u '1 2]
+ (ax dxdy


fJ J 02dxdy
0o o0

Once again using the solvability condition we obtain

--e = h + h + (I2
2 0 1 2 22

2 v(kG -G s)(1+B)tanha s
a op sc 5
2 21)Se 0 2
(k+B) (na +1/2vtanha s)
s s



32 vG
h 2 -(3 2a +I )A 2- I Se v
1 a 22 23 21 0 2

(kG c-G S)(a + v/2tanha )/tanha.


I2 v(G c- ) k(a + v/2 tanha )
2 2 (k+B)2 tanhaz

(na -v/2 tanha s)B
2 So+ 121
n tanha s

Se v(kG9 -G /n)
(k so (B).3-11)

and B (a -1/2 vtanha )tanha s/(na + 1/2 vtanha s)tanha
S& s s s s.

L2 3

2 ax2
o ax



4.4. Calculations and Comparisons

Using the results of the previous section, Se2 was calculated for

various values of the experimental parameters and the results are shown

in Figs. 4-6 to 4-9. Figure 4-6 was drawn in order to check the

derivation with that of Ungar and Brown (1984a) and shows calculations

done for 2-D rolls with Ds 0 and 8 = 1 but is a more complete

calculation for the experimental parameter ranges involved than

theirs. The calculations agree with those of Ungar and Brown but it is

an unexpected result nevertheless, showing multiple regions of

subcritical and supercritical instability. It also seemed to contradict

the earlier calculations of Wollkind and Segel (1970) who did not see

any supercritical instability but this was resolved in the paper of

Wollkind and Wang (1988) and hence Fig. 4-6 agrees with their

calculations as well. As argued in Section 3.4 if we treat the amin

curve as the experimental "operating line," the bifurcation is mostly

supercritical which also seems to contradict the experimental results of

de Cheveigne et al. (1986) who observed only subcritical bifurcation.

In Fig. 4-7 calculations were done for more realistic values of Ds

and 8 and the results are significantly different with only one region

of supercritically and another of subcriticality. Ignoring solid

diffusion introduces a singularity to the problem and this accounts for

the distortions observed in Fig. 4-6. In Fig. 4-7 if we move along the

operating line for a fixed liquid temperature gradient, initially for

small growth velocities the bifurcation is supercritical, until at a

critical velocity a transition point is reached and the bifurcation

becomes subcritical. Hence for every imposed liquid temperature




4000 -


Figure 4-6.

50 100 150 200 250 300
VELOCITY (gms'1)

Velocity vs. wavenumber chart for the one-sided
approximations of Ungar and Brown for 2-D rolls


20 40 60 80

VELOCITY (gms-1)

Figure 4-7.

Velocity vs. wavenumber chart for 2-D rolls
using the Pb-Sn system

gradient there will be a critical growth velocity below which the

bifurcation for 2-D roll disturbances is always supercritical. The

really surprising result here is that when these transition points for

various values of the gradient are joined, we get a straight line

through the origin. We now have a clearly demarcated supercritical

upper-triangular and subcritical lower-triangular one. It is well known

that the onset condition Seo does not change much in the neighborhood of

amin (see for example Coriell, McFadden and Sekerka (1985)) and so amin

is more of an interval than a unique point. It can also be seen from

Fig. 4-7 that the amin curve practically hugs the line of transition

from sub to supercriticality and if we impose an interval for amin it

would straddle the transition line. Thus for 2-D rolls it is unlikely

that we will ever see a sharp transition from subcritical to

supercritical bifurcation in experiments. More likely we will observe

subcritical behavior throughout as reported by de Cheveigne et al.

Proceeding to the three dimensional patterns, we obtained almost

identical results for square cells and cross rolls. Figure 4-8 is for

square cells and we see quite a change with the supercritical region

acquiring a characteristic balloon shape and having a sharp transition

to subcritical bifurcation along the amin curve. But here too if these

points of transition are joined, a straight line is the result,

demarcating a supercritical upper-triangular operating region and a

subcritical lower triangular one. Unlike the case of 2-D rolls these

should be visible to the experimentalist. So in order to avoid cross

rolls and square cellular instabilities not only should the crystal be

grown when Se < Seomin, but one should do so in the upper triangular

region. Figure 4-9 demonstrates the universality of our result in being

0 2 4 6 8 10 12 14
VELOCITY (gms-1)

Figure 4-8.

Velocity vs. wavenumber chart for square cells
using the Pb-Sn system

-60 -40 -20

VELOCITY (gms-1)

Figure 4-9.

Velocity vs. wavenumber chart for the fusion
problem for square cells using C-Austenite system






applicable for the melting problem as well and repeating the derivation

for this case separately as was done by Wollkind and Raissi (1974) is


Finally, even though for hexagons Se1 was usually non-zero (and

hence Se2 cannot easily be calculated), in Section 4.2 we saw that there

were points at which Se1 did go to zero. If we attempt to evaluate Se2

for hexagons at these points I21, 122 and 123 become

21 15 (1 + p2) (4.4-1)

5 2 2
122 a (1+p) (4.4-2)
22 4

23 3 a (1+p2) (4.4-3)

and as can be expected Se2 will be in the form

Se = Se (1 + p2) (4.4-4)

where Se2 is that corresponding to p = 0. The other possible value is

when p is /3 or -/3. Here Se2p can be calculated from (4.3-9) (4.3-

12) and (4.4-1) (4.4-3). As mentioned in Section 4.2 these points are

usually far away from amin and in Fig. 4-5 we found that along this

curve Se2 is positive at low values of "a." As "a" increases, at a

point above a=amin, Se2 becomes negative. Hence along this line the

bifurcation diagram is forward bending at low "a's" (including amin) and

becomes backward bending at high values of "a." Even though this is

true only along the Se1=0 line, by analogy with other cell patterns we

conjecture that this is valid when Se *0 as well. In other words we

expect the bifurcation diagram to be forward bending along the ain line

and below (as shown in Figs. 4-10 and 4-11) but a transition to backward

bending along the amin line could occur at high velocities. Far above

the amin line the bifurcation diagram should be backward bending (see

Figs. 4-12 and 4-13).

To confirm these conjectures we turned to the nonlinear

calculations for hexagons of McFadden et al. (1987) but unfortunately

they were unable to complete the bifurcation diagram as their attempts

to compute the curve for <0O failed. Also their calculations were done

only for the case of p=0 and not for p=/3. But they did confirm the

existence of subcritical instability.

To sum up, it has been shown that the Mullins and Sekerka and the

Woodruff models of morphological instability are of limited validity.

The uniform formulation is the more exact representation of the problem


it is applicable for all growth velocities and not just the

relatively rapid solidification and fusion regions and provides a

single formulation from which all the different models for

various growth conditions can be obtained as limiting cases,

eliminating the duplication of derivations for different cases;

it incorporates the whole crystal and melt regions into the

problem and not just a boundary layer region adjoining the

interface facilitating the study of various domain effects like

convection on morphological instability;

-- Unstable

Sei =
Se2 >


, SeSO


P =

^ -p =+

Figure 4-10.

Bifurcation diagram
for hexagonal cells

Figure 4-11.

Bifurcation diagram
for hexagonal cells



- Stable
-- Unstable

Se e
/ Sei = 0
Se2 > 0O

p =1= 0

Figure 4-12. Bifurcation diagram
for hexagonal cells

Se, 0

SP = 0

Figure 4-13. Bifurcation diagram
for hexagonal cells

it avoids the incorrect predictions of subcritical bifurcation

regions because of the singularities inherent in previous models.

The principle of exchange of stabilities has been shown to be

applicable to this model as well even though only in an

asymptotic sense.

When the weakly nonlinear technique of Malkus and Veronis (1959) is

applied to this problem in a systematic way, it resulted in important

information about the shape of the bifurcation diagram for various

growth conditions. Some of these results are similar to those obtained

for Marangoni instability (see Joseph (1976)) which leads us to assert

that these results are valid for all hydrodynamic instability problems

in which the nonlinearity lies only on the boundary.

If rectangular cells, cross rolls or two dimensional rolls are

the horizontal cell patterns then the bifurcation diagram will

always be locally symmetric. For hexagonal cells or cylindrical

rolls they are generally non-symmetric and hence the bifurcation

is subcritical.

The problem can display a multiplicity of solutions for the same

eigenvalue. Specifically for hexagons there are two possible


Considering the morphological instability problem in particular the

following were shown:

for two-dimensional rolls there are two operating regions, one

subcritical (backward bending) and the other supercritical

(forward bending), but since the demarcation is not sharp its

probable that only subscritical bifurcation will be observed



* For rectangular cells the forward bending region has a

characteristic balloon-like shape and here too there is a

straight line dividing the operating region into subcritical and

supercritical zones, but here the transition is sharp and hence

probably observable experimentally.

* For hexagonal cells and cylindrical rolls the bifurcation diagram

shows both backward and forward bending behavior but the exact

regions of each can only be conjectured.


5.1. Rayleigh-Marangoni Convection in Brief

When a horizontal layer of quiescent fluid is heated from below, on

account of thermal expansion, the fluid at the bottom will be lighter

than the fluid at the top. This unstable arrangement is maintained by

the viscosity of the fluid which inhibits any flow and suppresses

disturbances such that there will be a stable conduction profile in the

fluid. But when the adverse temperature gradient exceeds a critical

value, the viscous force is overcome and there will be cellular

convection. This instability is known as Rayleigh-Benard convective


There are several variations of this problem. Instead of an ad-

verse temperature gradient, there could be an adverse solute concentra-

tion gradient in the fluid causing once again an unstable top-heavy

arrangement. The convective instability arising from this is known as

solutal convection or the solutal Rayleigh problem. Another way to

cause convection is to have a very thin fluid layer heated from below,

but the top surface of the fluid instead of being kept at a fixed lower

temperature, is allowed to remain free. Here the thinness of the fluid

layer makes buoyancy effects negligible but convection will be caused by

surface tension variation on the free surface. This is known as Maran-

goni convection. Finally, there could be combinations of the above

three types of convective instability. When convection is caused by

thermal and concentration gradients it is known as double-diffusive

convection. Combination of either thermal or solutal convection with

surface tension driven flow is the Rayleigh-Marangoni problem.

In addition to these there are several other combinations possible

like Rayleigh-Benard convection with rotation or with a magnetic field

but for purposes of comparison with morphological instability it will be

seen that the three causes for convection mentioned above are the most


In the discussion to follow it is desirable to consider the most

general form of this problem. Despite there being several interesting

features in the problem of double-diffusive convection, the causes for

convection there, the temperature gradient and the solute concentration

gradient are both bery similar and it is sufficient to look at the

effect of one gradient. The manner in which the surface tension effects

the problem is very different from the buoyancy effects and a general

formulation should include both. Accordingly we will examine the Ray-

leigh-Marangoni problem with thermal convection. In the following

sections several other reasons for looking at this problem will become


The equations for the Rayleigh-Marangoni problem are given by Sarma

(1987) and we will reproduce them here and refer the interested reader

to 'his paper for details. The steady state dimensionless Boussinesq

equations in the domain are

V*V = 0


V p + a Tg V*VV = 0

V2T PrV*VT 0



where p is the modified pressure, g the acceleration due to gravity, Ra

the Rayleigh number and Pr the Prandtl number. The boundary conditions

at the bottom of the fluid layer are

at z = 0,

T = T and w = 0, aw/3z = 0


where w is the vertical component of velocity.

The boundary conditions at the top are more complicated. Not only

will there be surface tension variation across the surface but the

surface is also free to deflect like the liquid-solid interface in

crystal growth.

at z = 1 + C,

VT-n = BiT

V-n = 0

- Bopn + Crn*[VV + VVT] = MaV HT + fin


The dimensionless quantities are Bi the Biot number, Bo the Bond number,

Cr the Crispation number, Ma the Marangoni number and H the surface

curvature (see Scriven and Sternling (1964) for details).

Initially there will be a quiescent, linear, stable, conducting

solution to the problem with V = 0. At a critical value of the charac-



teristic parameter (Ma or Ra) this conducting solution becomes unstable

to infinitesimal perturbations and we have a convective solution.

Performing a linear stability analysis about the conduction state,

separating variables and doing considerable manipulations we get in the


(D a2 38 = a Ra6

2 2 3 2
(D a ) w = a Raw



where "a" is the wavenumber, w

component of velocity and 6 the

the Fourier coefficient of the vertical

Fourier coefficient of the temperature.

At the boundary at z = 0,

w = Dw = 8 = 0

At z = 1,

w = 0




BiD w = a MaD8

BiCr(D w 3a2Dw) = a2(Bo + a2)(De + Bi9)


When the density variation with temperature is negligible or in the

absence of gravity, then Ra = 0 and we have the pure Marangoni problem

with all the nonlinearities only in the boundary. Even with this effect

in the boundary the important result of bifurcation, namely the fluid

velocity, effects only the domain, the deflections in the boundary being

only a secondary effect of convection. On the other hand, when the

upper boundary is also kept at a fixed temperature we have, Ma = Bo = Cr

= 0 and so

at z = 1, w Dw = e 0 (5.1-14)

This then is the pure Rayleigh-Benard problem with all the nonlin-

earities only in the domain and the resulting nonquiescent solution also

manifests itself in the domain as convection. The Rayleigh-Marangoni

problem described by eqns. (5.1-8) (5.1-13) is a mixed problem with

nonlinearities in the domain and the boundary but the convective solu-

tion resulting from these nonlinearities shows up mainly in the domain.

5.2. The Augmented Morphological Problem

As can be seen from Section 5.1, in the Rayleigh-Marangoni problem

there is a Ra-Ma domain-boundary duality which does not seem to exist in

morphological instability. From the problem description in Chapter 3 it

is easy to see that all the nonlinearities for this problem lie only in

the liquid-solid interface. This is a limitation because by virtue of

being on the boundary the Sekerka number is unique and hence also has a

unique eigenfunction and is insufficient when solutions to inhomogeneous

versions of the linearized morphological instability problem are needed,

as in Chapter 6 where "imperfections" are considered. This difficulty

also crops up in the pure Marangoni problem, but the Rayleigh-Marangoni

problem comes to the rescue, as there are countably many corresponding


values of Ra for each value of Ma and hence also countably many

eigenfunctions forming a complete set (see Rosenblat, Homsy and Davis

(1981)). The naturally occurring duality of Ma and Ra enables solutions

to inhomogeneous problems to be obtained in a straightforward manner.

In morphological instability there is no such obvious, naturally

occurring boundary-domain duality and it is necessary to create one. To

avoid confusion we will refer to the pure morphological problem of

Chapter 3 as the Sekerka problem, and (by analogy with the Rayleigh-

Marangoni problem) set up an eigenvalue problem, with the eigenvalue in

the domain, which we will call the "augmented morphological problem."

liquid: V2p = MpL (5.2-1)

2 H
S- v 3q- = (5.2-2)

solid: Vp2 = Mp (5.2-3)

2 s
nVq 3 = 0 (5.2-4)

where M is the eigenvalue which we label as the morphological number.

The boundary conditions are

at z = -1, pt qt = 0 (5.2-5)

at z = s, p = 0 (5.2-6)

pa pS + t(G Gs) =0 (


at z = 0,

pt + Seqe + t(G + SeG t) + A (5.2-8)

I 8 a 0 (5.2-9)
az az

kqg qs + t(kG o G ) = 0 (5.2-10)

n 3- vq + vq = 0 (5.2-11)
Tz 'az Z

with periodic conditions at the lateral boundary at r = R. We can now

separate variables expressing the horizontal dependence as zeroth order

Bessel's functions of the first kind Jo(air/R), where a, are the zeros

of the first order Bessel's functions of the first kind J1.

pt(r,z) = E Z P (z) J (a.r/R) (5.2-12)
i=1 j=1 0 1

In addition if we take Fourier transforms in the horizontal direc-

tion and solve for q.., q .. and t the equations reduce to a system

in pij and p ij. If we define a column vector

Qij = (5.2-13)

and a matrix differential operator Li

2 2
(D a )
L. = (5.2-14)
2 2
BY.(D a.)
1 1

then the domain equations reduce to

1 ij -ij 1Qj

where Y =(G + Se aA )/(G + Se G a2 )
i o i s c i

G = (G B G)/(B + k)

(a v/2tanha )tanha .s
-1 9_i s31
S (na + v/2tanha .s)tanha z
Si 31 Ui

2 2 /4 1/2
ai =(a.i + v /4)

2 2 2 1/2
a = (a2 + v /4n )

The boundary conditions become

at z = -1,

at z = s,

at z = 0,


p = 0


= 0

BiQij = 0

Bi D











Defining an inner product

0 -0 a -
P tij P2ik + jPsij Psik dz
-1 0


where the "*" refers to the adjoint eigenfunction and "-" the complex


It can easily be seen that the system described by eqns. (5.2-13) -

(5.2-24) is self-adjoint in this inner product and so the eigenfunctions

Qij are complete. Solving the system we get


A Sin X7 a7 (z + 1)
Aij ij i

A Sin/M. (s z)
sij ij i


where Mij are solutions of the equation

Se = (- GT + aA )/ Go
T 1


k G tan/M a S
T k tanM a2S +
Z ij i

+ k G tan/M .- ai

k G tan/M. a-
ss ij i

Here A ,. and Asij can be determined from the normalizing condition
xtij sij


= 6jk

where 6 is the Kronecker delta.


5.3. Comparison of Morphological Instability with
Rayleigh-Marangoni Convection

The principal aim of this section is to relate the mathematical

characteristics of the two problems so that we may introduce some of the

extensive mathematical techniques used to study Rayleigh-Marangoni

convection to morphological instability, but we will make some physical

comparisons as well. The augmented morphological problem described in

Section 4.2 is similar to Rayleigh-Marangoni convection. The augmented

problem is self adjoint but the Rayleigh-Marangoni problem is nonself

adjoint. Both have an infinity of eigenvalues and corresponding eigen-

functions, but while completeness of the eigenspace is assured for the

former, special theorems are required to show this for the latter (cf.

Nadarajah and Narayanan (1987)). It should also be noted that while the

Rayleigh-Marangoni problem attempts to describe a realistic situation,

the augmented morphological problem was artificially created in order to

solve inhomogeneous versions of the Sekerka problem described in Sec-

tions 3.1 and 3.2.

This brings us to the question whether there is a practical situa-

tion which is described by this mathematical concoction. The difficulty

in coming up with one stems from another important difference between

the two problems. In the pure Rayleigh-Benard problem (where Ma, Bo and

Cr are all zero) the nonlinearity is in the domain and the instability

too manifests in the domain as convection. Even in the pure Marangoni

problem (where Ra is zero) where the nonlinearity is in the boundary,

the instability is still mainly in the domain. In contrast, in the

Sekerka problem the nonlinearities and the resulting instability show up

in the boundary. Though there are other boundary effects (like kinetic

undercooling) which can cause morphological instability the only domain

effect which could give M physical significance is a heat source term in

the form MT or Me-E/RT (see Joseph (1965)). We do not know of any

experiment where morphological instability was observed as a result of a

heat source in the melt or the crystal but if one does exist it will

provide the true analogy to Rayleigh-Benard convection.

This is relevant as Hurle (1985) has attempted a comparison between

the Rayleigh-Benard problem and the Sekerka problem. It can now be seen

that the Sekerka problem can only be compared to the pure Marangoni

problem, with Se corresponding to Ma and A corresponding to the reci-

procal of Bo. Besides, for periodic lateral boundary conditions, the

eigenvalues of both problems, Se and Ma, are unique.

Based on this comparison we can make an important conjecture.

Vrentas, Narayanan and Agrawal (1981) have shown that for the Marangoni

and the Rayleigh-Marangoni problems when the nonperiodic no-slip condi-

tion for velocity is imposed at the sidewalls, the eigenvalue Ma is no

longer unique and has countably many values. In other words when the

walls are a finite distance apart Ma has many values, but as they are

gradually moved apart we have "spectral crowding" and in the limit when

they are sufficiently far apart to impose the periodic boundary condi-

tion of total slip, all the values of Ma coalesce into a unique num-

ber. Recently, following Coriell et al. (1980), several workers have

looked at the coupled problem of morphological instability with solutal

convection and all have assumed periodic boundary conditions. We sus-

pect that here too if the no-slip condition for velocity at the side

wall is imposed, the Sekerka number will no longer be unique. All this

raises the question of completeness of the Marangoni and the Sekerka

eigenspaces and its probable that generalized eigensolutions (see Nai-

mark (1967)) are needed when Ma and Se are chosen as eigenvalues.

When these two problems are considered in a finite container we can

see yet another difference. Both problems have simple eigenvalues

except at certain aspect ratios of the container where two horizontal

modes can coexist (cf. Rosenblat, Homsy and Davis (1982)). In a typical

experiment of Rayleigh-Benard convection we would expect to see a dozen

or so convection cells (see for example Koschmieder (1967)) and increas-

ing or decreasing the number of cells by one can significantly effect

the problem. Hence the multiple points in this problem are extremely

important and have been the subject of study. But in morphological

instability a single alloy crystal can contain hundreds of individual

cells and the addition or loss of one has hardly a noticeable effect on

the problem and consequently multiplicity of the lateral eigenfunctions

loses its significance. In addition unlike the Rayleigh number it is

well known that near the critical wave number aN, the critical value of

the Sekerka number Seo hardly changes (see for example Coriell, McFadden

and Sekerka (1985)) and the choice of aN has very little effect on

SeoN. Conversely, the choice of the operating Se will have a tremendous

impact on the resulting wavenumber (cf. Ramprasad and Brown (1987)).

Other differences have been mentioned in Chapter 4. Both the

Marangoni problem and the Sekerka problem have symmetric bifurcation

diagrams near the bifurcation point for two-dimensional rolls and rec-

tangular cells and nonsymmetric curves for hexagonal cells and cylindri-

cal rolls. But in morphological instability the curves can be "forward

bending" or "backward bending" depending on the operating conditions,

whereas in Marangoni convection the curves are forward bending every-


where. Hence the occurrence of subcritical instability is more wide-

spread in morphological instability.


6.1. Nature of Imperfections

When the morphological instability problem was formulated in Chap-

ter 3, several effects were ignored and the resulting problem is an

idealized or "perfect" one. Inclusion of these can alter the problem in

several ways, for example kinetic undercooling of the melt becomes an

important effect in rapid solidification, but all it does is alter the

onset condition for morphological instability. In the parlance of

bifurcation theory, an "imperfection" is an effect on the "perfect"

problem which alters it in a specific way. Such an imperfection will

cause the morphological instability problem not to have a planar solu-

tion at all even below the onset condition, and these are known as

bifurcation breaking imperfections.

The effect of a typical imperfection on the bifurcation diagram is

shown in Fig. 6-1. The broken line is the solution in the presence of

imperfection and it can be seen that the interface will be nonplanar for

all nonzero values of Se. Obtaining solutions to the problem with

imperfections is extremely difficult and we will only seek asymptotic

solutions. Hence the problems to be considered should have very small

imperfections. Under such conditions the method of matched asymptotic

expansions of Matkowsky and Reiss (1977) can be employed and here it

will be used in a way similar to the work of Tavantzis, Reiss and Mat-

kowsky (1978) for the Rayleigh-Benard problem.
-? I

- -

Se =Seo N

Figure 6-1.

Imperfect bifurcation diagram showing inner
and outer expansions

The fairly straightforward. The variables are expanded

asymptotically with the imperfection parameter about the planar and the

nonplanar solutions and two outer expansions 00O and 01 are obtained as

shown in Fig. 6-1. At the bifurcation point SeoN these expansions break

down and it is necessary to have inner expansions I1 and 12 near SeoN

and to join the corresponding 00 and 01, matching conditions have to be


The imperfections that can be analyzed in this fashion must, of

course, be small effects else a full-blown nonlinear solution will be

needed. Two of the most important effects which are habitually ignored,

namely imperfect insulation of the ampoule wall and advection in the

melt, are such imperfections and readily lend themselves to this type of

analysis. In Chapter 3 when the morphological instability problem was

modelled by a uniform formulation, it was mentioned that one reason for

this was to consider a finite crystal/melt region. This finiteness was

only in the vertical direction and in the horizontal direction the

imposition of periodic boundary conditions effectively meant that the

ampoule side walls were infinitely far apart. The two imperfections

that are to be considered are caused by nonperiodic lateral boundary

conditions and another way of looking at the effect of these

imperfections is to say that the container is now being considered to be

finite in the lateral as well as the vertical direction.

6.2. Imperfection Due to Heat Loss

As mentioned in the last section, it has been customary in this

problem to assume periodic boundary conditions laterally, which is


equivalent to assuming that the walls of the ampoule are perfectly

insulated or that they are so far apart that their effect can be ig-

nored. In practice neither of these is likely to be achieved and here

we will examine the effect of a small amount of heat loss or heat gain

from the wall on the problem. If we take the ampoule to be cylindrical

with radius R,

r = R, liquid: -k2 6fD (z) (6.2-1)

solid: -k --s = 6f (z) (6.2-2)

where f and f. are such that f (S) = fs () and fX(-1) = f (s) = 0.

If f and fs are positive it will mean heat loss and if they are nega-

tive, heat gain. If we make the transformations

T (z,r) = 6 k r + 9 (z,r) (6.2-3)

f (z)
and T (z,r) = 6 k r + 9 (z,r) (6.2-4)
s k s

and substitute these in the steady-state versions of eqns. (3.1-15) -

(3.1-24), the temperature equations in the domain become

V2 = 6rD2f /k (6.2-5)

V2 = 6rD2f /k (6.2-6)
s 5 s

The outer boundary conditions will remain unchanged but at the interface

eqn. (3.1-21) becomes

6- = 6r(f /k f /k ) (6.2-7)

6 + SeC + AH = 6rf /k (6.2-8)

and eqn. (3.1-22) converts to

7V Cn 8V6 *n = Lv dr(Df Df ) (6.2-9)

In order to solve this system we will be treating the heat loss as

an imperfection on the perfect problem (i.e. when 6 = 0). The perfect

problem is of course the Sekerka problem of Section 3.1.

6.3. The Outer Expansions

As this problem has been defined in a finite geometry, the number

of cells are fixed by the container size and the growth conditions. We

will choose these so that aN is very close to amin the wavelength corre-

sponding to Seomin the least value of SeoN. Another important decision

is the selection of the wave pattern and our analysis is done for cylin-

drical rolls. An objection to this could be raised on the grounds that

in most experiments it is the hexagonal pattern which is observed. We

justify our assumption by once again making a comparison with Rayleigh-

Benard convection. It appears that in Rayleigh-Benard convection the

wave pattern selection is strongly influenced by the container size and

shape, with the hexagonal pattern prevailing for all shapes in wide

containers while in narrower ones the container shape determining the

pattern (e.g., cylindrical rolls for circular containers. See Kosch-

mieder (1967)). Since the principal aim of this paper is to study a

finite geometry effect on morphological instability, the cylindrical

roll pattern would be the logical choice. This is especially valid for

experiments such as those of Peteves (1986) where ampoules of radius

0.025 inches were used.

In Chapter 4 it was shown that the bifurcation diagram is unsymme-

tric for cylindrical rolls and the form of the outer expansions 00 and

01 are shown schematically in Fig. 6-1. If we use superscript o to

identify the problem with perfect insulation we can seek solutions by

means of asymptotic expansions about the perfect problem.

k = 6 (6.3-1)

with similar expansions for s, CE, Cs and ;. Substituting these expan-

sions into the problem and collecting the terms of order 6 we get an

inhomogeneous linear system. If we separate variables horizontally

1 1
9= 6 J (a r/R) (6.3-2)
t 1t o i

i=1 1

and eliminate C., C and C we get for the expansion about the planar
jL s

L 1 f= (6.3-3)
1 i i




1 2 I

i f /k
s s

. = 2j rJ (a.r/R)dr/ rJ2(a.r/R)dr
0 0 1

The boundary conditions are

at z = -1,

at z = s,

at z = 0,

a =0

6 =

1 1
1ii 1

h = I
f1 i /

-Y.f /k
1 s s


The eigenvalue problem of Section 5.2 has been shown to have a

complete set of eigenfunctions and hence can be used to solve the above











91 E= Q .J (a.r/R) (6.3-11)
i=1 j=1

J(Q h) ..,.
ij i M M.
ij z=0 1J

where J(Q .,h )|1 is the bilinear concomitant evaluated at z=0
ij i z=0

J(Q = P- (Df Df )Ii/k, + DPi (f/k Y f /k )I.
ij 0iis s i


1 1
If we assume that hN and f are nonzero, then when Se = SeoN, MN1

becomes zero and the outer expansion (6.3-11) (6.3-13) will fail.

When we expand about the nonplanar solution we end up with a much

more complicated system. But since we have expressed the solution to

the perfect nonplanar problem itself in terms of a perturbation series

in E, we can do the same for the imperfect problem.

1 1 1
= + Eo + ... (6.3-14)
o 1

If we take the zeroth order problem and once again separate vari-

ables laterally,

1 1
S E J (a.r/R) (6.3-15)

1 01 0 1

and eliminate of Se Ctio and 5.2, e get the same system described by

eigenfunctions of Section 5.2, the solution is

1 1
0o = E Q J (a.r/R)
i=1 j=1


with taking the same form as eqns. (6.3-12) and (6.3-13).

Since these equations are valid only for Se close to SeoN, we can write

an expansion for MN1

M (Se) = M (Se ) + (Se Se) dM
N1 Ni oN oN dSe =SeON

Lmplifies to

MN1(Se) = 0 + eSelNMN1 + 0( 2)

N dSe Se = Se

+ ... (6.3-17)



and SelN has been derived in Section 4.2. Hence the outer expansion

about the nonplanar solution is

9 = + *- J (a r/R) + 0(e2) +
c oN o N

( J)
6(Q N> -- J (a r/R) + 0(E )] + 0(62)

which si



6.4. The Inner Expansions

The outer expansions described in the

SeoN, requiring inner expansions in this


Se(y) = Se + Slb +

6(u) = Ce

previous section fail at Se =

vicinity. For Se very near

k kb/k!
E /k!



where b and c are integers to be determined such that the solution does

not fail at Se = SeoN. We can now expand the solutions in orders of p.

k k
u(r,z,p) = 6(r,z,Se,6) = Z u (r,z)ip /k!

v(r,z,p) = C(r,z,Se,6) = E v k(r,z) k/k!

w(r,z,i) = S(r,Se,6) = Z w (r)k /k!




Here the superscripts on u, v and w are not powers but indices.

Substituting these expansions into the problem set up in Section 6.2 and

collecting the terms of order 1j, we once again get an- inhomogeneous

problem. Separating variables and eliminating vQ, vs and w, the problem

reduces to

L. 1 = 6'f1 (6.4-6)

with boundary conditions

at z = -1,

at z = s,

at z = 0,


u 0

U = 0

Bi = 6'h

9 -

6, d S dSe
d6 =0 Se d =0


The variables Li and Bi have been defined in Section 5.2 and f.
a 1
and h. in Section 6.3. But Yi now becomes
1 1

2 2
Y. = (G, + Se c G a2A)/(G + Se G a2 A)
i z oN c i s oN c i


This means that the operator Li is singular when i=N, j=1 and hence

for solutions to exist 6'5f has to orthogonal to the null space of LN.

Using this solvability condition with QN1 we get

6'[J(QN1, N) z=0

- ^N1'N> = 0






From this we can conclude that 6' = 0 (assuming J ) and

that c>1. This makes the system (6.4-6) (6.4-12) homogeneous and

so are constant multiples of 9 .

9 = At


Proceeding to the next order we have

2 1
L.i = 6"f

and at z=0,

2 2
B.. = h.
11 1

G 1o Gc/(G + Se G -ca.A
h2 = 6"h1 2ASe' o1 c z oN c i + 2A
1 1 0 )1


S= rJ3(a r/R)dr/ rJ2 (a.r/R)dr
0 0





and Q is a complex function of the linear stability variables.

solvability condition for i=N, j=1 now gives


J(Qh) 6 1 = 0
Nl'hi) 'QN1c'fNr>=

Hence we can take 6" # 0 and Se' 0 and so b=1 and c=2. Equation

(6.4-19) is a quadratic in A with two real roots A, and A2 of opposite

signs, corresponding to two inner expansions Ii and 12.

Se = Se + 6 1/2 + 0(6)

= o + A 61/2 + 0(6)

In terms of 1, the

eqn. (6.3-11) can now be

outer expansion 00 about the planar state of

expressed by

1 = ( A JO(a r/R)/M^
*N1 N1 N N1 0o(aNr/R)/gM 1

and the nonplanar expansion 01 of eqn. (6.3-20)

N1 Nl'> J)/MNT + D/SelN]QN1 J(ar/R)
IO1 [


where we expressed oN as DNQN1 J (aNr/R). In attempting to match these

with the inner expansions of eqn. (6.4-21) we get

lim A =lim A2 1QN' f > J)/SM1

lim A = lim A2 = D S/Se 1
1 +2 N4 -N



Hence using Al in eqn. (6.4-21) will give the inner expansion I1

matching 00 for Se < SeoN with 01 for Se > SeoN. A2 will result in

matching 00 for Se > SeoN with 01 for Se < SeoN as shown in Fig. 6-






6.5. Imperfection Due to Advection in the Melt

In modelling morphological instability during solidification it is

customary to neglect the difference in density between liquid and

solid. Though this difference is very small (e.g. for the Pb-Sn system,

the ratio of densities p /pa is 1.035) the volume contraction upon

solidification results in closed streamlined flow in the melt which

fundamentally alters the planar state of this problem. Obtaining the

solution to the planar state with this flow is difficult enough even

without attempting the formidable task of solving the nonplanar prob-

lem. In order to include the effect of the flow on morphological insta-

bility previous workers have therefore relied on simplifying assumptions

or looked at limiting cases.

The scenario in a typical solidification experiment is schemati-

cally shown in Fig. 6-2. The heating coils maintain constant tempera-

tures at positions z = -Z and z = s in the ampoule. As solidification

proceeds, the ampoule .is pulled in the positive z-direction at a velo-

city v in order that the liquid-solid interface will remain stationary

at z = 0. To fill the space created by volume contraction upon solidi-

fication, the melt will move towards the crystal with a bulk velocity

of v(p /p- 1) and this process is commonly referred to as "advec-

tion." In the presence of the rigid ampoule walls, this flow also

resembles the bolus flow of slugs through a pipe, resulting in closed

streamlined flow in the melt.

In the literature we find three approaches to tackling this prob-

lem, the least of which is the work of Caroli, Caroli, Misbah and Roulet

(1985b) who ignored the rigid side walls. Then the only effect of

- Encapsulant

--- z=-Im

z= -




Crystal growth with advection in the melt

Figure 6-2.

advection is a negligible modification of the growth velocity in the

melt. The other two approaches are more substantial and analize

limiting cases of this phenomenon. Since the traditional formulation of

morphological instability concentrates only at a "boundary layer" region

of the liquid-solid interface, in such a model the closed streamlined

flow shown in Fig. 6-2 can be approximated by Couette flow in Region I

and stagnation flow in Region II. Delves (1968, 1971 and 1974) and

Coriell, McFadden, Boisvert and Sekerka (1984) looked at the effect of a

forced plane Couette flow and showed that the onset of morphological

instability is somewhat suppressed for disturbances in the flow

direction, while disturbances perpendicular to the flow were

uneffected. Recently McFadden, Coriell and Alexander (1988) examined

the effect of a plane stagnation flow on disturbances perpendicular to

the flow and here too the flow was found to be stabilizing.

Based on these two limiting cases it might be tempting to conclude

that advection generally stabilizes the liquid-solid interface. In this

and the next few sections we will embark upon a more complete analysis

of advection than has been attempted so afar and show that the above

assertion is questionable at best. In our model we consider morphologi-

cal instability with advection in the absence of natural convection. We

choose not to include natural convection as we wish to study the effect

of advection only and natural convection will only further complicate an

already nontrivial problem. Besides it has been shown very conclusively

by Coriell, Cordes, Boettinger and Sekerka (1980) and by Caroli, Caroli,

Misbah and Roulet (1985a) that the convective and morphological modes

are decoupled except at points where they become comparable. Hence our

model will be valid for the high growth velocity region where morpholo-

gical instability is dominant and also in low gravity environments where

natural convection can be neglected.

The experimental set up was briefly described earlier. The steady

state domain equations in dimensionless form are

V2T = 0


nV C v as = 0
3 3z





Here u is the velocity of the closed streamlined flow in the melt

caused by advection and 6 is (p S/p 1). The boundary conditions are

at z = -1,

at z = s,

at z = ;,

TQ = T1, C = C1

TS = T2

T = T = SeC +
2. s 2

VT *n BVT *n = Lv
X. 3

kC = C

VC,*n (1+6)vCe = nVCs*n vC
X> JO S.







Before we can write the equations for u, further assumptions are

necessary. The domain of the velocity equations extend over the entire

melt region, not just the depth X that we have considered. In the float

zone technique this melt region will still have a constant depth as

solidification proceeds but in Bridgman growth the melt region will

continuously shrink. Since crystal growth velocities are usually so

small, even in the latter technique it takes awhile before there is an

appreciable change in the melt depth. Hence over a short time span a

constant depth assumption will be valid even for Bridgman growth. We

will also assume that the melt is bounded on all sides by rigid walls.

This will require a container for zone refining and a very viscous

encapsulant for the melt in Bridgman growth.

Under these assumptions the velocity problem has been solved by

Duda and Vrentas (1971) for low velocities and here we will only give

the solution and refer the interested reader to their paper for the


u = (6vR/r)3I/3z, u = 6v(1 (R/r)3i/ar) (6.5-11)
r z

S= Z A F (r)Sinnr(1 + z/Y) + Z G (z)rJ (b r/R) (6.5-12)
n n n=1 n 1 n
n=1 n=1

F (r) = [r (nirr/Y)I2 (n/h) r21 (n7/h)I (nrr/Y)/R /12(n/h) (6.5-13)

G n(z) = 2Bn[Sinhb n(z+Y)/R (z/Y+1)exp(-b nz/R)Sinhb nh]
n n n n n

+ 2C (z+Y)/R expb h Sinhb z/R
n n n


9 = J2(b )[4lb2hexpb h + (2expb h exp(-b h) exp3b h)/h]
n o n n n n n n

Bg = -4b hexpb h Z A Q +
nn n n m mn

2(1 exp2b h) E A Q (-1)m'
n mn

hg C = 4(b hexpb h Sinhb h) Z A Q -

(4b h+2 2exp2b h) E A Q (-1)m

bn n-r[I (nw/h)I2(nw/h) I 2(nr/h)]/I2(nr/h)

2b (mit/h)2 2(mw/h)J2(b )
mn (mh)[(/h2 b2
mn = I (mit/h)[(mw/h)2 + b T
2 n

Y = am /Z,






h = Y/R

where Im is the total melt depth, R the dimensionless radius of the

ampoule, ur and uz the radial and axial components of the velocity, Jp

the Bessel function of the first kind of order p and Ip the modified

Bessel function of the first kind of order p. Duda and Vrentas have

determined the first twenty coefficients of An for various values of h

and we will use these in our calculations.


6.6 Nonexistence of the Planar State

The key assumption made in the two previous approaches to solving

the advection problem was that a stable planar state was still possible

even in the presence of melt flows. But this is true only if the velo-

city profile very close to the interface is assumed to be in a special

form, the vertical component of velocity being independent of any hori-

zontal coordinates and the horizontal component of velocity at most a

linear function of that horizontal coordinate. As can be seen by eqns.

(6.5-11) (6.5-19) these assumptions can only be viewed as limiting

cases and in general ur and uz will be complex functions of both r and

z. In this section we will use this to prove the nonexistence of planar

solutions for this problem. This proof can also be used to show the

nonexistence of planar solutions for the case of imperfect insulation

considered in Sections 6.6 6.4.

In Section 3.1 where we assumed that 6 was zero, a planar solution

was possible with the temperatures and concentrations functions of z

only. Then the energy balance condition at the interface (given by eqn.

(6.5-8)) could be used to determine the constant growth velocity v. In

other words the choice of T1, T2 and C1 fixed v. For the case of non-

zero 6 we will show that planar solutions do not exist by contradic-

tion. If there is a possible planar solution eqn. (6.5-2) specifies

that C cannot have radially independent solutions. The solution will

be in the form

C (r,z) = E C (z) J (a r/R) (6.6-1)