Citation

Material Information

Title:
A study of free and moving boundary problems involving thin crystal growth
Creator:
Liang, Shin-Jye, 1958-
Publication Date:
Language:
English
Physical Description:
xvii, 193 leaves : ill. ; 29 cm.

Subjects

Subjects / Keywords:
Bond number ( jstor )
Boundary conditions ( jstor )
Convection ( jstor )
Crystal growth ( jstor )
Crystals ( jstor )
Geometric angles ( jstor )
Heat transfer ( jstor )
Hysteresis ( jstor )
Simulations ( jstor )
Solidification ( jstor )
Aerospace Engineering, Mechanics, and Engineering Science thesis Ph. D ( lcsh )
Boundary value problems ( lcsh )
Dissertations, Academic -- Aerospace Engineering, Mechanics, and Engineering Science -- UF ( lcsh )
Genre:
bibliography ( marcgt )
non-fiction ( marcgt )

Notes

Thesis:
Thesis (Ph. D.)--University of Florida, 1993.
Bibliography:
Includes bibliographical references (leaves 183-192).
General Note:
Typescript.
General Note:
Vita.
Statement of Responsibility:
by Shin-Jye Liang.

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
Resource Identifier:
030010770 ( ALEPH )
29809577 ( OCLC )

Full Text

A STUDY OF FREE AND MOVING BOUNDARY PROBLEMS
INVOLVING THIN CRYSTAL GROWTH

By

Shin-Jye Liang

A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF
THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY

UNIVERSITY OF FLORIDA

1993

To my parents, my wife, Susan, and my son, Kevin

ACKNOWLEDGEMENTS

I owe an immense debt of gratitude to my parents for the support they have

always offered and the love that never wavered. Their belief in me has always been a motivating force behind my achievements. And my wife, Susan, has helped me in more ways than she knows. providing constant support and a few lessons on what is and what should be important in life.

I thank my advisor, Professor Wei Shyy, for introducing me to this interesting field of research and providing me with his guidance. His patience and insights have greatly contributed to the quality of this work.

Thanks are due to the other members of my supervisory committee, Professors Bruce F. Carroll, Chen-Chi Hsu, Ulrich H. Kurzweg and Craig Saltiel, for their comments and encouragement throughout my time at the University of Florida.

Some experimental information used in this study was supplied by

Dr. D. G. Backman and Dr. D. Y. Wei; their help is greatly appreciated. I wish also to thank all my colleagues in the CFD group with whom I had helpful discussions, in particular, Dr. Renwei Mei and H. S. Udaykumar for many fruitful suggestions.

The present work is partially supported by GE Aircraft Engines, by the National Center for Supercomputer Applications at the University of Illinois, and by the National Science Foundation.

iii

Page

ACKNOWLEDGEMENTS

LIST OF TABLES ............................................

LIST OF FIGURES ...........................................

NOM ENCLATURE ..........................................

ABSTRA CT ................................................

CHAPTERS

1 INTRODUCTION .......................................

. . . . iii

. vi . vii xiii . xvi

1.1 Background ............
1.2 Description of Edge-Defined 1.3 Scope of the Present Work
1.4 Outline ...............

Film-Fed Growth Technique . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . .

2 INTERACTION OF TIME STEPPING AND CONVECTION

2.1 Introduction ....................
2.2 Model Problem Analysis ...........
2.3 Results and Discussions ...........
2.3.1 Linear Burgers Equation ....
2.3.2 Nonlinear Burgers Equation 2.4 Von Neumann and Spectral Analyses
2.5. Summary ......................

SCHEMES ....... 13

....... 13
....... 16
....... 22
....... 23
....... 25
....... 30
....... 33

3 MENISCUS CHARACTERISTICS WITH APPLICATION TO THE EFG
PROCESS .............................................. 53

3.1 Introduction .......................................... 53
3.2 Configuration, Assumptions and Formulation ................. .56

iv

. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .
. . . . . . . . . . . .

3.3 Energy Form ulation .................................... 60
3.4 N um erics ............................................ 62
3.5 Results and Discussion .................................. 62
3.5.1 Basic Meniscus Behavior with Upward Pulling .......... 63
3.5.2 Effect of Pull Direction on Meniscus Behavior .......... 67
3.5.3 Quasi-equilibrium Motion with Hysteresis ............. .70
3.6 Conclusions .......................................... 77

4 NUMERICAL TECHNIQUES FOR SOLVING FREE AND MOVING
BOUNDARY PROBLEMS WITH APPLICATION TO THE EFG
PRO C ESS .............................................. 95

4.1 Introduction .......................................... 95
4.2 Literature Review ..................................... 97
4.2.1 Exact Solutions ................................. 97
4.2.2 Approximation Solutions .......................... 97
4.2.3 Numerical Solutions ............................. .98
4.3 Front Tracking and Moving Grid ......................... .103
4.3.1 Assessment of the Quasi-stationary Approximation ...... 104 4.3.2 A General Procedure for Interface Tracking .......... 109
4.3.3 Assessment of Time Stepping Schemes for Phase Change
Problem s ................................... 111
4.4 EFG Process Physics .................................. 113
4.5 M athematical M odel .................................. 115
4.6 Numerical Simulation of the EFG Process .................. 124
4.7 Conclusions ......................................... 130

5 DYNAMIC SIMULATION OF THE EFG PROCESS ............ 146

5.1 Introduction ......................................... 146
5.2 Dynamic Perturbation on the EFG Process with St =1 ......... .148
5.2.1 Base Temperature Perturbation .................... 148
5.2.2Pull Speed Perturbation .......................... 151
5.3 Dynamic Perturbation on the EFG Process with St = 0.024 ..... .155 5.4 Dynamic Contact Angle ................................ 157
5.5 Summary and Concluding Remarks ........................ 160

6 CONCLUSIONS AND RECOMMENDATIONS ................ 178

6.1 Achievements and Findings ............................. .178
6.2 Directions for Future Research ........................... 181

REFERENCES .............................................. 183

BIOGRAPHICAL SKETCH ..................................... 193

v

LIST OF TABLES

Tables Page

2-1 Amplification Factors of Backward Euler Scheme for Linear
Burgers Equation ................................... 35

2-2 Amplification Factors of Crank-Nicolson Scheme for Linear
Burgers Equation ................................... .36

4-1 Assessment of time stepping scheme for phase change problem.
(a) Comparisons of the predicted interface location at r = 0.2 using backward Euler and Crank-Nicolson time stepping schemes for Eqs.
(4.3) and (4.4) with 11 and 21 grids.
(b) Comparisons of the required CPU time on MicroVax 3100 using backward Euler and Crank-Nicolson time stepping schemes for Eqs.
(4.3) and (4.4) with 11 and 21 grids ..................... .132

4-2 Thermophysical properties and physical dimensions of sapphire
(A1203) used in the present EFG simulations (Case 1: St = 1 and
Case 2: St = 0.024) ................................. 133

4-3 Definition and magnitude of key dimensionless parameters based on
values given in Table 4-1 (Case 1: St = 1, Case 2: St = 0.024)
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134

4-4 Two different scaling procedure and resulting nondimensional
energy equations and boundary conditions ................ 135

vi

LIST OF FIGURES

Figures Page

1-1 Schematic illustration of the EFG sapphire fibre growth equipment
used to produce 25 fibre with diameters in the range 0.025 - 0.15 mm. The inset shows a magnified view of the die tip and meniscus
region where fibre solidification occurs ................... .11

1-2 Schematic illustration of the heat transfer mechanisms in the EFG
process ........................................... 12

2-1 Comparison of solutions of unsteady linear Burgers equation using
different convection schemes, backward Euler time stepping,
P = 10, C = 1, and 21 nodal points ...................... 37

2-2 Comparison of solutions of unsteady linear Burgers equation
obtained using the backward Euler and Crank-Nicolson time
stepping schemes .................................... 38

2-3 Comparison of the backward Euler and Crank-Nicolson time
stepping schemes with four convection schemes for unsteady linear
Burgers equation using 81 nodal points, and P = 10. C = 1 . ... 39

2-4 Effect of time stepping schemes on solution accuracy of the
nonlinear Burgers equation with v = 5x103 , Ax = 0.02, and
a t = 0.01 ......................................... 40

2-5 The L-2 error norm with respect to Ax of the nonlinear Burgers
equation with v = 5x103 and At/Ax = 0.5 ................ 41

2-6 Effect of time stepping schemes on solution accuracy of the
nonlinear Burgers equation with v = 5x105 , Ax = 0.02, and
A t = 0.01 ......................................... 42

2-7 The L-2 error norm with respect to Ax of the nonlinear Burgers
equation with k = 5x105 and At/Ax = 0.5 ................ 43

vii

2-8 Effect of time stepping schemes on solution accuracy of the
nonlinear Burgers equation with v = 5x10- , Ax = 0.02, and
At = 0.05 ......................................... 44

2-9 Performance characteristics of the combined backward Euler and
first-order upwind schemes with periodic boundary condition,
P = 103 , C = 0.5,and Ax = 0.05 ....................... 45

2-10 Performance characteristics of the combined backward Euler and
second-order central difference schemes with periodic boundary
condition, P = 103 , C = 0.5,and Ax = 0.05 ............... .46

2-11 Performance characteristics of the combined backward Euler and
second-order upwind schemes with periodic boundary condition,
P = 103 , C = 0.5,and Ax = 0.05 ....................... 47

2-12 Performance characteristics of the combined backward Euler and
second-order QUICK schemes with periodic boundary condition,
P = 103 , C = 0.5,and Ax = 0.05 ....................... 48

2-13 Performance characteristics of the combined Crank-Nicolson and
first-order upwind schemes with periodic boundary condition,
P = 103 , C = 0.5,and Ax = 0.05 ....................... 49

2-14 Performance characteristics of the combined Crank-Nicolson and
second-order central difference schemes with periodic boundary
condition, P = 103, C = 0.5, and Ax = 0.05 ................ 50

2-15 Performance characteristics of the combined Crank-Nicolson and
second-order upwind schemes with periodic boundary condition,
P = 103 , C = 0.5,and Ax = 0.05 ....................... 51

2-16 Performance characteristics of the combined Crank-Nicolson and
second-order QUICK schemes with periodic boundary condition,
P = 103 , C = 0.5,and Ax = 0.05 ....................... 52

3-1 Schematic of meniscus corresponding to edge-defined film-fed
growth process ..................................... 79

3-2 Characteristics of meniscus formation for Bo = 4x l04 , AP = 0,
h,/rb = 0.5, and upward pulling ......................... 80
(a) variation of 4. with angle Ob .
(b) permissible range of r/rb with respect to O.
(c) profile shapes for different angles Ob .

viii

(d) profiles corresponding to 0, = 60* .
(e) E profile v.s. 0, , for the specified value of 4O = 600
(f) profiles corresponding to 4, = 900 .
(g) E profile v.s.0. , for the specified value of I. = 90*.

3-3 Characteristics of meniscus formation for Bo=2.5xl04' , AP=5%,
hc/rb = 2, and upward pulling .......................... 82
(a) profiles corresponding to 00 = 16* .
(b) E profile v.s.4 , for the specified value of 4 = 16* .
(c) profiles corresponding to 0, = 20' .
(d) E profile v.s. 0, , for the specified value of 4, = 20' .

3-4 Characteristics of meniscus formation for Bo=2.5x103 , AP=5%,
h,/rb = 2, and upward pulling .......................... 83
(a) profiles corresponding to 0. = 60' .
(b) E profile v.s. 4e , for the specified value of 4 600 .
(c) profiles corresponding to 4 = 90' .
(d) E profile v.s. 0e , for the specified value of 4, = 900 .

3-5 Characteristics of meniscus formation for Bo=4x104 , AP=5%,
he/rb = 5, and upward pulling .......................... 84
(a) profiles corresponding to 4, = 90' .
(b) E profile v.s. , , for the specified value of 4, = 90* .
(c) profiles corresponding to 4, = 1300 .
(d) E profile v.s. , , for the specified value of 0, = 1300 .

3-6 Effect of the direction of pulling on the meniscus characteristics
with Bo=102 and AP = 0 ............................. 85
(a) permissible range of contact angle.
(b) permissible range of trijunction location.
(c) permissible meniscus shapes at 0, = 100 .
(Note that with given Bond number and hc/rb less than 1, the
meniscus shape is dominated by the geometrical constraint from the boundary conditions. Hence the meniscus shapes are essentially the same between two pulling directions. The effect of pulling direction
becomes noticeable only when hC/rb becomes larger.)

3-7 Effect of pressurization on the meniscus shapes with Bo= 102 and
AP = 1%. (excess interior pressure) ..................... .86
(a) permissible range of contact angle.
(b) permissible range of trijunction location.
(c) permissible meniscus shapes at 0(, = 700 .
(Note that the hydrostatic pressure effect becomes noticeable only
for a higher hc/rb , as evidenced by difference caused by pulling
directions.)

ix

3-8 Energy profiles for Bo = 10-2 , AP = I%, , = 700 and various
he/rb. Both upward and downward pulling cases are shown . ... 87

3-9 Schematic presentation of hysteresis model ................ 88

3-10 Meniscus behavior at lower Bond number, with no hysteresis.
( Bo = 10- , AP = 5%. h,' = 05 = 100) .......... 89
(a) Meniscus profiles obtained for the duration of oscillation.
(b) Variation of height h. with time.
(c) Computed value of top radius rc versus time.

3-11 Meniscus behavior at lower Bond number, with small hysteresis
range. ( Bo = 10- , AP = 5%. hc' = 0.5, 6A = 950, 6, = 85*) .. 90
(a) Meniscus profiles obtained for the duration of oscillation.
(b) Variation of height h, with time.
(c) Computed value of top radius rc versus time.
(d) Computed variation of contact angle 4 versus time.

3-12 Meniscus behavior at higher Bond number, with no hysteresis.
(Bo = 2.5x10' , AP = 5%. h, = 0.5, 6, = 6R =_100"1 ....... 91
(a) Meniscus profiles obtained for the duration of oscillation.
(b) Variation of height h. with time.
(c) Computed value of top radius r. versus time.

3-13 Meniscus behavior at higher Bond number, with no hysteresis.
( Bo = 2.5x10l.,AP = 5%. h = 1.5,A = 6- = 100 ....... .92
(a) Meniscus profiles obtained for the duration of oscillation.
(b) Variation of height h, with time.
(c) Computed value of top radius rc versus time.

3-14 Meniscus behavior at higher Bond number, with hysteresis.
(Bo = 2.5x10 ', AP = 5%. h = 1.5,d)A = 1100, 6R = 100) ... 93
(a) Meniscus profiles obtained for the duration of oscillation.
(b) Variation of height h, with time.
(c) Computed value of top radius r. versus time.
(d) Computed variation of contact angle d, versus time.

3-15 Meniscus behavior at higher Bond number, with small hysteresis
range. ( Bo = 2.5x10' , AP = 5%. hV = 0.5,dA, = 900, b, = 880) 94
(a) Meniscus profiles obtained for the duration of oscillation.
(b) Variation of height hc with time.
(c) Computed value of top radius r, versus time.
(d) Computed variation of contact angle 0c versus time.

x

4-1 Exact solution for melting of a solid confined to a semi-infinite
region ........................................... 137

4-2 Comparison of interface trajectories predicted by the coupled
(full equations) and decoupled (quasi-stationary equations)
approaches with three values of the Stefan numbers ......... .138

4-3 Comparison of computed and exact solutions for St = 0.1303,
St = 1.2216 and St = 2.8576. The figures on the left are the exact
and computed temperature fields at different time instants. The
figures on the right show the superposed exact and computed
interface locations versus time ......................... .139

4-4 Schematic configurations and boundary conditions of the present
EFG process ...................................... 140

4-5 Comparison of linear meniscus profile with the Laplace-Young
solution with RI = 1, H, = 0.188, R, = 0.763, and various values of Bond number. (a) Ap = 5 % and (b) Ap = 0 % ........... 141

4-6 The grid system and isothermal contours of the initial conditions for
both Stefan numbers ................................ 142

4-7 The grid system and isothermal contours of the steady-state solution
for Case 1: St = I .................................. 143

4-8 The grid system and isothermal contours of the steady-state solution
for Case 2: St = 0.024 ............................... 144

4-9 Effect of scaling procedure on the convergence characteristics for
Case 2: St = 0.024, choice 1: VN = O(St) and choice 2:
VN = (1) ....................................... 145

5-1 Time histories of H, and R, subject to a single harmonic
perturbation of Ob with (a) Q = 500 Ar and (b) 0 = 20 AT
(St = 1.0) ........................................ 163

5-2 Time histories of Hc and R, subject to a three-harmonic
perturbation of 0. with Q, = 500 Ar, '2 = 292 Ar, and
03 = 108 Ar (St = 1.0) .............................. . 164

5-3 The interface shapes at different time instants of (a) a singleharmonic perturbation . in 6b with Q, = 500 Ar and (b) a threeharmonic perturbation in ob with 0, = 500 Ar, 02 = 292 Ar, and
03 = 108 Ar (St = 1.0) .............................. 165

xi

5-4 Characteristics of interface movement with respect to different
forcing frequencies of Ob (a) The relative phase angles between
HC and 0,, and between H, and R. . (b) Percentage variations of
H, and R, (St = 1.0) ................................ 166

5-5 Time histories H, and R_ subject to a single harmonic perturbation
of UP with (a) Q = 500 Ar and (b) ( = 20 Ar (St = 1.0) . ... 167

5-6 Characteristics of interface movement with respect to different
forcing frequencies of U, (a) The relative phase angles between He and Up, and between H, and R.. (b) Percentage variations of H.
and R, (St = 1.0) .................................. 168

5-7 Time histories He and Rc subject to a three-harmonic perturbation
of U, with Q, = 500 AT, (2, = 292 AT and (2 = 108 Ar (St = 1.059

5-8 Time histories of saw-toothed forcing with Q = 20 Ar and
A = 0.1, and the corresponding H, and R, (St = 1.0) ....... .170
(a) base temperature perturbation Ob(T),
(b) pull speed perturbation U,(T).

5-9 Time histories of H. and R, subject to a single harmonic
perturbation of U, with (a) period=500 AT and (b) period=20 Ar
(St = 0.024) ...................................... 171

5-10 Characteristics of interface movement with respect to different
forcing periods of U, (St = 0.024) ...................... 172

5-11 A typical Bode diagram depicting sensitivity of the fibre diameter
variation with respect to different forcing frequencies of U,
(St = 0.024) ...................................... 173

5-12 Model of dynamic contact condition at the trijunction point
(St = 0.024) ...................................... 174

5-13 Time histories of H. and R, subject to a single harmonic
perturbation of U, with period = L00Ar (St = 0.024) ........ 175

5-14 Time histories of dH/dr and dRlIdr subject to a single
harmonic perturbation of U, with period = LOOAT (St = 0.024) 176

5-15 Characteristics of interface movement with respect to different
forcing periods of U, for different static contact angles
(St = 0.024) ...................................... 177

xii

NOMENCLATURE

Ai
a A-B/ B.E.
Bi Bo
C C.-N.
Cp
E E, E2 E3
e2 EFG
F

f
f(x,4) G, G,
g HC
h

h,
i
i,j Km
k
L Lm
1

Ma
N

integers denoting grid point in space coordinate wavenumber of the m-th component thermal conductivity nondimensional 1 wavelength of the m-th component dimensional vertical length of the physical domain latent heat of fusion reference length scale Marangoni number nondimensional unit normal vector

amplitude of external forcing of the i* harmonics
capillary length
Backward Euler time stepping scheme
Biot number
Bond number
Courant number
Crank-Nicolson time stepping scheme
specific heat at constant pressure
total free energy
gravitational potential energy
surface energy of meniscus
surface energy of the wetted area of the top plate
L-2 error norm
process Edge-Defined Film-fed Growth process
1) nondimensional f(z,t)
2) dimensionless function defining interface x = s(y,t)
meniscus shape distribution as function of r and t
right hand side of Eq. (2.16)
amplification magnitude
heat flux along the top of the solid phase [
gravitational acceleration
nondimensional h,
1) height of the melt/solid interface h(r,t)
2) convective heat transfer coefficient [W/
height of the trijunction point
V-I

[cm-']
[W/cm-K]

[cm] [cm] [J/g] [cm]

xiii

[cm]

[J/g-K]
[erg] [erg] [erg] [erg]

[cm]

W/cm2] [cm/s2]

[cm]
cm2-K]
[cm]

n
Nr Nz
p

Pe qj, q2, q3, J
R
r
Ra Rad rb r, rf (r,z,t) (R,Z,t')
S
s
St
T
t
TA Tb
To Td
Tm t,
TVD UP
(U,),vg
u

VN
W , W,
x
(x,y,t) (X,Y,t)
z
z

unit normal vector along the boundary total number of grid points along the r- direction total number of grid points along the z- direction 1) cell Peclet number
2) pressure [dyne/cm2]
Peclet number
metrics coefficients defined in Eq. (4-44) nondimensional radius
Rayleigh number
radius of the trijunction point [cm]
roughness factor
physical coordinates
nondimensional physical coordinates nondimensional interface shape S(R,Z,r) = Z - H(R,r) = 0 dimensional function defining interface [cm]
Stefan number
temperature [K]
temporal coordinate [s]
ambient temperature [K]
base temperature [K]
truncation error for approximating convective term truncation error for approximating diffusive term melting temperature [K]
reference time scale [s]
total variation diminishing nondimensional pull speed nondimensional average pull speed convective velocity [cm/s]
dimensional pull speed [cm/s]
normal velocity of interface [cm/s]
nondimensional velocity components in z- and r- direction [cm/s] spatial coordinate [cm]
physical coordinates
nondimensional physical coordinates vertical coordinate [cm]
nondimensional vertical coordinate

Greek

a

thermal diffusivity thermal expansion coefficient

[cm2/s]
[1/K]

xiv

F
-Y

Ap Ap
E
6

(Ob)avg

P

OA
sR Ob

(OPr)
p

nondimensional thermal diffusivity surface tension [dyne
infinite small increment density difference between liquid and gas phases [g/
pressure difference between liquid and gas phases [dyne/ emissivity
1) phase angle
2) nondimensional temperature average base temperature root of Eq. (4-10) dynamic viscosity [g/
kinematic viscosity [c
transformed coordinates density [g
Stefan-Boltzman constant [5.670x 10-12 W/cm
1) unknown in Burgers equation 2) instantaneous contact angle of the meniscus at the trijunction
point
advancing contact angle receding contact angle shooting angle at the die tip contact angle
contact angle at steady-state apparent contact angle period of external forcing of the i* harmonics

i-cm] m2/s]

/cm3] 2-K4]

[deg] [deg] [deg] [deg] [deg] [deg] [deg]
[s]

nondimensional quantity

ambient quantity gas phase liquid phase melting point previous time step current time step next time step radius component solid phase vertical component

xv

/cm]

/cm3] cm21

[deg]

superscript

*

subscripts

a
g
1
m
n-1
n n+l
r
s
z

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

A STUDY OF FREE AND MOVING BOUNDARY PROBLEMS INVOLVING THIN CRYSTAL GROWTH by

Shin-Jye Liang

May 1993

Chairperson: Wei Shyy
Major Department: Aerospace Engineering, Mechanics, and Engineering Science

Thin sapphire fibers, which find applications in high temperature composites and optical sensors, can be grown from the melt by the Edge-defined Film-fed Growth (EFG) technique. In this thesis, a thermocapillary model is developed to simulate and illustrate the crystal solidification characteristics. The model, which is axisymmetric and appropriately scaled, solves the energy equation in both the solid and liquid phases separately, and tracks the movement and evolution of their interface explicitly using a newly developed method based on a combined Lagrangian/Eulerian scheme.

The formation of meniscus profiles is examined utilizing the concept of minimization of free energy. Meniscus behavior influenced by Bond number,

xvi

pressurization, aspect ratio. contact angle, and pulling direction has been investigated using both static and quasi-dynamic models for the contact condition.

In the course of dynamic simulation, the characteristics of the EFG system

subject to base temperature perturbation and pull speed perturbation are studied. To both types of perturbation, the solid/melt interface responds at the forcing frequency. Due to the appearance of multiple time scales in the solidification process, both the model prediction and the experimental observation indicate that as the forcing frequency increases, its effects on crystal growth are reduced. Furthermore, the interface behaviors yielded by the static and dynamic models for the contact conditions at the trijunction point are different, indicating the inadequacy of the conventional approach. The capability developed and insight gained from this work can help improve the EFG process.

xvii

CHAPTER 1
INTRODUCTION

1.1 Background

Free and moving boundaries are encountered in a variety of engineering applications. Analysis of such problems is often difficult because the free and moving boundaries, whose locations are unknown a priori and must be determined as a part of the solution, make the problem nonlinear. The interaction between phase change and other thermofluid transport processes adds additional complexities, making such problems intrinsically difficult to solve. For example, in crystal growth from the melt, the coupled heat transfer and fluid dynamic equations, conjugated with phase change, solid phase heat conduction and radiative heat exchange with the ambient surfaces, must be solved simultaneously in order to understand the characteristics of the process, the behavior of the growing crystal and the effects of other manufacturing process parameters that determine the quality of the crystal (Tseng et al. 1992, Brown 1988, Viskanta 1988).

Mathematically, the free and moving boundary problems are challenging and, except for a few simple cases, have to be solved numerically. Considerable progress has been made in the solution of problems involving free and moving surfaces using numerical methods., A central problem in the numerical modeling

1

2

of such processes is the treatment of heat and mass flow conditions that arise at the moving solid/liquid interface and the liquid/gas meniscus whose shapes and locations are unknown and vary in time. Many approximation techniques and numerical methods have been developed, as reviewed by Floryan and Rasmussen (1989), Crank (1984), Ozisik (1980), Wilson et al. (1978), Bankoff (1964), Carslaw and Jaeger (1959).

The characteristics of the crystal/melt interface is dictated by the Stefan number since it determines the growth rate of the crystal into the melt. In order to control the shape/quality of the crystal, the combined effects of fluid (melt) motion and heat transfer in a variable domain bounded by capillary (melt/gas) and Stefan-type (solid/melt) interfaces need to be studied. Open research problems are of both mathematical (existence, uniqueness of solutions) and practical natures. The main purpose of this study is to numerically analyze free and moving boundary problems involving phase change with special emphasis on the edge-defined film-fed growth (EFG) process for producing thin sapphire crystals.

1.2 Description of Edge-defined Film-fed Growth Technique

Sapphire fibers for use in both optical sensors and structural composites have been produced using the edge-defined film-fed growth (EFG) process (Kalejs et al. 1983ab, LaBelle 1980, Sachs 1980, Swartz et al. 1975), wherein a crystallographically oriented fibre is grown from a meniscus of molten alumina.

3

An inert die is used to control the extent of the melt and to help shape the melt and thus the resulting crystal. The EFG process can yield continuous lengths of sapphire fibre with diameters in the range of 0.05 to 0.25 mm. The fibre grower, schematically shown in Figure 1-1 (Backman 1992), contains three main components: the hot zone, the puller system, and the fibre spooling system. The hot zone consists of an inductively heated refractory metal crucible, with the edgedefining dies all enclosed within a water-cooled, environmentally regulated chamber. The melt is supplied to each die tip from a capillary-fed manifold. The puller system consists of a belt puller and a fibre guide. The belt puller utilizes a double belt traction mechanism driven by a precision stepper motor and is monitored with an optical encoder. After the fibre exits the puller, it passes through several pulleys and is wound on spools under regulated tension. Within the hot zone, the fibre growth process begins with seeding at the die tip, which establishes a crystallographic orientation. The puller speed is generally kept constant; the meniscus dimensions and fibre diameter are manipulated through changes in the induction coil power level setting to adjust die tip temperature.

In practice, both the die tip temperature and the puller speed may

experience either intentional or unintentional variations in time. The unintentional variations can result from the perturbations to the operating conditions caused by fluctuations in the environment, coil power fluctuations, or motor speed irregularities. The intentional variations can be designed to compensate some of the aforementioned perturbations to achieve a controlled growth process and uniform crystal properties and dimensions.

4

Throughout the development of the EFG process, several problems have been encountered that have potentially severe effects on the quality of the grown crystal. For example, the shape/quality of the thin crystals grown from the melt is strongly influenced by the intricate coupling of heat transfer and melt flow in the EFG system. The heat transport environment must yield a well-controlled temperature field in both the solid and liquid phases so that the crystallization rate, the shape of the solidification interface, and thermoelastic stresses in the crystal can be controlled. A general goal of the design of heat transfer systems for any crystal growth configuration is to establish a nearly constant temperature gradient along the crystal/melt interface and to control the cooling rate of the crystal. Low dislocation and defect densities occur when the temperature gradients in the crystal are low (Kalejs et al. 1983a,b, Sachs 1980 and Swartz et al. 1975). However, because of the small dimensions of the thin crystal sheets, which are of the order of a few mm, the rate of heat losses to the surroundings along the surface of the crystal/melt could be large. This usually results in large temperature gradients near the solidification boundary, which are undesirable. To overcome this problem, the shape-defining dies can be enclosed within a watercooled, environmentally regulated chamber, as illustrated in Figure 1-1. Active control may be desirable to produce crystals of uniform cross section because the shape of the crystal is constrained only by capillary action.

A schematic illustration of the heat transfer mechanisms in the EFG

process is shown in Figure 1-2. A film of melt forms from the top of a die and is

5

solidified by convection and radiative heat transfer to the surroundings. It is important to control the growth process according to the heat transfer characteristics in the system to obtain a crystal of definite diameter and desired homogeneity. In the melt, convection driven by mechanisms such as pulling of the seed crystal, buoyancy and surface tension will affect the heat transfer process and hence the solidified crystals. For thin crystal growth from the melt via the EFG technique with a dimension of a few mm, the order of magnitude of convective heat transfer is expected to be relatively small compared with that of the conductive one. Order of magnitude analysis of these heat transfer mechanisms under real operating conditions are evaluated in Chapter 4. Based on this analysis, a mathematical model is developed for simulating the dynamic characteristics of the crystal growth process. Furthermore, the issue of an appropriate scaling procedure, which plays a critical role in determining the computational efficiency, will be addressed.

1.3 Scope of the Present Work

The purpose of this study is to mathematically analyze the transport

phenomena in the EFG system so as to establish the relationships between the operating conditions and the crystal shape/quality. To develop a complete capability of numerical simulation of free and moving boundary problems involving phase change, several technical issues need to be addressed, including the assessment of time stepping, convection and diffusion schemes, meniscus

6

behavior and the contact condition at the trijunction location, the solution algorithm for the thermal field in both liquid phase and solid phase, the front tracking technique, and the adaptive gridding to accommodate time-dependent and irregular domain configurations. In order to develop a comprehensive model to conduct dynamic simulations of the EFG process, the following aspects are considered in this thesis.

First, the discretization schemes for the temporal and spatial terms and

their interaction in the time-dependent thermofluid computations are investigated. In order to clearly understand the issue, both linear and nonlinear Burgers' equations are used as model problems to study the numerical performance of the various schemes. The von Neumann stability and fast Fourier transform techniques are utilized to analyze the stability, accuracy, and dispersive and diffusive characteristics of the combined performance of time stepping and convection schemes. The scope of this aspect of the study is more generic and is intended to provide some general guidelines for helping devise suitable computational strategies for growing crystals of large sizes, where convection effects are more pronounced and may vary spatially.

Second, the physics and appropriate treatment of the solid-liquid-gas

trijunction point, which is not well understood at the present time, is explored. Consideration of force equilibrium leads to a stress singularity at the contact line. The modeling of static contact lines, where a contact angle is specified at the trijunction point, is justified only in the liquid-gas-solid system under the

7

assumption that the solid wall is smooth and homogeneous. The extension of this concept to dynamic contact lines is desirable, to relate the instantaneous contact angle to the velocity of contact line motion. The attempts to determine this function from physical experiments have not been very successful and this remains an open area of research. Relatively detailed discussion will be given later, for both the static and dynamic contact conditions. A simplified dynamic model, taking into account the effect of the direction and speed of contact point movement, will also be proposed and incorporated into the present model. Based on both the static and dynamic contact conditions, meniscus characteristics and solidification behavior in the EFG process will be examined and compared.

The third aspect of this study is to simultaneously solve the equations of

heat transfer in both solid and liquid phases, liquid/solid interface movement, and the dynamics of meniscus and contact lines using a curvilinear coordinate system. Body-fitted coordinates are used for mapping the irregular shape of the timewisevarying solid/liquid and liquid/gas interfaces. A new Lagrangian approach is developed for tracking the phase front, and a solution procedure is used for solving simultaneously the thermal fields in both liquid and solid phases. This yields a complete EFG numerical simulation.

Finally, thermal sensitivity and stability of the EFG for sapphire fibre

production is numerically studied. The dynamic behavior of the EFG system with respect to an external disturbance, e.g.,varying the base temperature or pull speed, is investigated. The proposed dynamic model is incorporated to study the

8

dynamic EFG characteristics. Based on this model, the variation of radius of the grown crystal is compared with the static wetting condition, which is conventionally employed in numerical simulations of crystal growth (Crowley 1983, Ettouney and Brown 1983). The relationship between response and forcing frequency, in terms of magnitude and phase shift, is identified.

1.4 Outline

The present study therefore can be divided into several parts:

discretization schemes for temporal and spatial terms and their interaction, meniscus and contact point behaviors, front tracking and solution techniques for phase change problems, and dynamic simulation of the EFG process.

In Chapter 2, the time stepping schemes, the convection schemes, and their interaction in unsteady thermofluid computations are presented. Both linear and nonlinear Burgers' equations are used as model problems to study the stability, accuracy, and dispersive and diffusive characteristics of each scheme. With the aid of a stability analysis, the overall performance of the scheme is discussed.

The shape of the melt/gas meniscus that connects the growing crystal with the die, which is governed by the Laplace-Young equation, is critical in determining the properties of the solidified crystal. In Chapter 3, the meniscus characteristics with application to the EFG process are studied. Both static and dynamic contact line motions are investigated. Quasi-equilibrium dynamics of the meniscus is simulated using a simplified hysteresis model for the contact angle at

9

the trijunction point. The implication of the variation of the meniscus shape on the crystal growth characteristics is explored in detail.

In Chapter 4, a general treatment of the interface movement is designed, which is not limited to a single-valued or isothermal surface employed by other investigators (Lacroix 1990, Thompson and Szekely 1989, Ettouney and Brown 1983, Crowley 1983, Sparrow et al. 1977). The development of the numerical

technique is reported in Shyy et al. (1992b), which demonstrates that the method is robust and can handle both isothermal and nonisothermal phase change problems. The technique is combined with a moving gridding procedure, determining the thermal field in both the solid and liquid phases and the evolution of both the solid/melt interface and the melt/gas meniscus to predict the solidification characteristics of sapphire fiber growth via the EFG technique for different Stefan numbers. Two different scaling procedures for nondimensionalizing the governing equations are examined, and the computational characteristics and relative efficiency of these two procedures are compared.

In Chapter 5, the sensitivity and dynamic characteristics of the EFG system subject to perturbations of the operating environment are explored. The dynamic behavior of the EFG system with respect to the external disturbance, i.e., variations in either base temperature at the die tip or pull speed, is investigated. A dynamic model is also incorporated into the system and compared with the conventional static model for the contact condition at the trijunction point. In Chapter 6, a summary and some suggestions for further work are given.

10

Overall, the contributions of this thesis are to develop and utilize advanced and original numerical techniques for simulating the EFG process, to examine and assess the contact conditions relevant to the present physical problem, and to construct a predictive tool incorporating both the physical models and the numerical techniques. The developed model system is envisaged as a computational tool, to make detailed investigations of the dynamic characteristics of the EFG processes.

G

Induction

0 to Spo Puller t

uide N

Coil a * le - - -

Fiber

nlpr

Crucible
Top Meniscus

Chamber Die Tip

View Port

Figure 1-1 Schematic illustration of the EFG sapphire fibre growth
equipment used to produce 25 fibre with diameters in the range
0.025 ~ 0.15 mm. The inset shows a magnified view of the die tip
and meniscus region where fibre solidification occurs.

{

Crucib

12

Figure 1-2

pulling speed

Solid
(diffusion)

interface
trijunction (latent heat) point

meniscus Liquid (Young-Laplace Eq)
diffusion & convection)

Die

Schematic illustration of heat transfer mechanisms in EFG process.

CHAPTER 2
INTERACTION OF TIME STEPPING AND CONVECTION SCHEMES FOR UNSTEADY FLOW SIMULATION

2.1 Introduction

The numerical solutions of heat transfer, fluid flow, and other related transport problems can be obtained by solving the laws governing these processes in discrete forms. By using appropriate numerical techniques to discretize the spatial derivatives, the original continuous problems are reduced to a system of ordinary differential equations in time. A time stepping scheme can then be used to advance the solution in time from given initial and boundary conditions. Categorically, two kinds of time stepping methods can be defined, namely, explicit and implicit. The explicit schemes provide a noniterative "marching" procedure for obtaining the solution of each nodal point at present time in terms of the known preceding and boundary values. Explicit schemes limit the time step size due to the stability as well as the accuracy requirements. Implicit procedures, on the other hand, generally involve the simultaneous solution of all unknowns in terms of both the known preceding and the unknown present values at other nodes. The implicit methods are usually more stable, and their restriction in time step size is mainly due to the accuracy consideration.

13

14

Besides the time stepping methods, the convection schemes play an important role in both steady and unsteady flow computations. Discretization of the convection terms is one of the principal difficulties in computing the complex fluid flows. The first-order upwind scheme used to be a popular choice due to its superior numerical stability as well as the wiggles-free characteristics. However, the resulting accuracy yielded by this scheme is often found to be inadequate in view of its first-order truncation error and strong numerical, as opposed to physical, diffusion effects (Shyy 1985, Roache 1972).

Many alternative convection schemes have been studied, as reviewed in Hirsch (1990), Fletcher (1988), Correa and Shyy (1987), and Warming et al. (1973), and the findings for the same schemes are not necessarily consistent among different studies; they are obviously problem as well as implementation dependent. Furthermore, it is known that the performance of a given convection scheme can be quite different between the steady and unsteady flow problems (Shyy 1984), hence, it is important to address the issues related to the interactions between the time stepping and the convection schemes for unsteady flow calculations. In this chapter, assessment of the relative merits of the various possible choices for solving the convection-dominated transport problems is made through a series of well defined cases.

The viewpoint taken here is that a numerical scheme preferably is of high order of accuracy in truncation error not only so that it can perform well for problems of smooth solution profiles, but, equally importantly, so that no over/undershoots or spurious oscillations that violate the conservation laws appear, in order to

15

be able to honor physical realizability. In the latter regard, the concept of total variation diminishing (TVD) has been found very helpful for constructing an appropriate convection scheme (Hirsch 1990). It is known that, with the use of the forward Euler scheme, the only convection scheme that can guarantee the TVD property across a sharp gradient or a discontinuity is first-order (Roe 1986). Here, through pragmatic testing and basic analysis, the intent is to gain insights of the performance of the schemes considered without insisting on the satisfaction of TVD. The reasons are twofold. No known technique is available to enforce TVD for a uniformly high order time stepping scheme; furthermore, for problems with substantial source terms, the concept of TVD does not apply (Hirsch 1990).

Three time stepping schemes, the first-order backward Euler, the second-order Crank-Nicolson, and the third-order Adams-Bashforth/Adams-Moulton (A-B/A-M) predictor-corrector schemes, are considered. In conjunction, four convection schemes, first-order upwind, second-order central differencing, second-order upwind, and QUICK (Leonard 1979), are examined using both the linear and nonlinear Burgers equations. The goal is to understand the dispersive and dissipative characteristics exhibited by each scheme and to offer guidance for utilizing these features to help improve the overall numerical accuracy. The Burgers equation, of the linear as well as nonlinear forms, is used as the test problem. Furthermore, the von Neumann stability analysis and the FFT spectral analysis will also be used to help understand the characteristics of each formulation.

16

The results to be presented in the following are a part of the attempt of developing the necessary models for predicting the crystal growth processes with different physical scales and operating conditions. As will become clear later, the convective effect is not dominant in the present EFG processes; however, the work conducted here will become essential for any system of substantial convective effects. Assessment of the time stepping schemes for phase change problem, which is critical for the problem under study, is presented in section 4-3.

2.2 Model Problem Analysis

The model problem considered is the unsteady one-dimensional Burgers equation

+ (u ) = (vl40), 0:5x: I,t >0, P =constant >0 (2.1)
at ax Ox ax

where if u is prescribed constant, the equation is linear, while if u is taken to be a function of 0, then it is nonlinear. Both cases are adopted here since they are relevant to practical applications under different circumstances. The boundary and initial conditions are given in the following

(I) Linear Problem ( u = 1 )

(i) boundary conditions

(0,t) = 1

(1,t) = 0

17

(ii) initial condition

O(xO) = for 0O!5x!50.2
{0 , othenvise

(II) Nonlinear Problem ( u = )

(i) boundary condition

4(Ot) = 0

t
4(1,t) =t
1 +( ) exp( ) to 4vt

(ii) initial condition

4(x,1) =
1 +( ) exp(-) to 4

where t = exp(1/8v).

The corresponding analytic solution (Lohar and Jain 1981, Benton and Platzman 1972) is
x

O(x,t) = t 2 (2.2)
1 +( ) exp( ) to 4vi In what follows, the second-order central differencing scheme is always adopted to approximate the diffusion term, i.e.,

18

2 4 2, + Oi- + T (2.3)

where Td , of the order of O[(Ax)2], is the truncation error.

As to the convection term, the following schemes are considered and tested:

(a) first-order upwind

(u4), - (u4),1 _ 24
8(uk1+ T for u > 0
a(uq5) _Ax C (2.4)
ax (u4) )j. - (uO) + T, for u < 0
Ax

where T,, of order of O[(Ax)], is the truncation error.

(b) second-order central differencing

a(u(0) _ (uk)1, - (u)i -1 + T (2.5)
ax 2 Ax

where T, is O[(Ax)2]

(c) second-order upwind

1 [ 3(uO), -4(uo),_ +(uo)i-2 ] + T , foru >0 a(uo) 2Ax (2.6)
ax -(UO)i.2+4(uo),., -3(u0), ] + T, for u < 0
2Ax

where T. is O[(Ax)2]

19

(d) QUICK

- [3(u)j,. +3(u) -7(ub), _, +(uO)i -2]+T,, foru >0
8(ub) = 8Ax (2.7)
ax [ -(u)i.2 +7(u4), . -3(uI), -3(u), _1] +T ,for u <0
8Ax

where T, is also O[(Ax)2].

When the first-order backward Euler method for the time derivative is used, the discretized form of Eq. (2. 1) can be expressed as follows:

(a) first-order upwind

(1+P) , - (2+P+f-) 4, + 0 = -- - 0 (2.8)
C C

where P = u Ax/v = cell Peclet number,

C = u At/Ax = Courant number

and the superscript "o" represents the known value.

(b) second-order central differencing

P P P P 0 29
(1+f) 4,_1 - (2+-) 4, + (1- ) 4,, = - - 4 (2.9)
2 C' 2 C'

(c) second-order upwind

- - + (1+2P) _ - (2+-P+.) 4, + 4,., = - 4 (2.10)
2 2 C C'

20

(d) QUICK

-P4, s-+(1 + 7 P)O, -(2+-P+)4,1+( 1- P)k +=-4,O (2.11)
8 8 8 c 8 C

The Crank-Nicolson scheme is a second-order accurate method and usually is described as unconditionally stable (Fletcher 1988, Anderson et al. 1984). However, as demonstrated later, this does not mean that the solutions are free of spurious oscillations. If the convection and diffusion terms are replaced by the appropriate discretization schemes at both the current and previous time levels, the corresponding discretized equations result:

(a) first-order upwind

(1+P) 4 - (2+P+21) P + 4k., = C (2.12)
- (1 +P) q5 + (2+P-2-) 4
C

(b) second-order central differencing

(1+P) 4,_ - (2+2)4 +(1- ) 4,) =
2 C 2 (2.13)
- (1+-) 0_ + (2-2.) P v - (l.-E) 0.
2 C ' 2 '+1

21

(c) second-order upwind

- Oi-2 +(I+2P) Oi - (2+3P+2-) O, + 4,=
2 2 C (2.14)

P (i- -+ (2+3P-2P) 4 - .
2 2 C

(d) QUICK

P-2 + (13+P) (2-+(2+P+2 ) 3 + ) Oi8 8 8 C 8 (2.15)
P 7 3 P 3 )0
Oi- -7 -(1 + P) 0- + (2 + P-2 ) 00 - 0- P) i+.
8 8 8 C 8

The third-order predictor-corrector scheme, based on the Adams-Bashforth scheme as the predictor step and the Adams-Moulton scheme as the corrector step, is explicit; it involves successive computations at each node but does not need to solve the whole domain simultaneously.

If Eq. (2.1) is rewritten in the following manner,

ao = - _ (u4) + a(P-4) =J(x,k) (2.16)
at ax ax ax

then the present predictor-corrector algorithm results Adams-Bashforth method :

0(00 = O(n) + A( 3f(") - f("-1 ) (2.17)

22

n'i) = (n) + A ( 5f(n+I* + 8f(n) - f(nl) ) (2.18) 12

where (n+1)* designates the value obtained from the predictor step. It is noted that the appropriate convection and diffusion schemes can be used to represent the function f(x,4k). The present A-B/A-M scheme needs a starting procedure since it requires information from the previous two nodes. Because the solutions close to the left boundary are smooth, it is found that no discernible difference exists between using the Crank-Nicolson and the backward Euler schemes as the starting procedure.

2.3 Results and Discussions

The viewpoint taken here is that a numerical scheme preferably is of high order of accuracy in truncation error not only so that it can perform well for problems of smooth solution profiles, but equally importantly, so that no over/undershoots or spurious oscillations that violate the conservation laws appear, in order to be able to honor physical realizability. In the latter regard, the concept of TVD has been found very helpful for constructing, as well as assuming, the convection scheme Hirsch (1990). It is known that, with the use of the forward Euler scheme, the only convection scheme that can guarantee the TVD property is first order. Here, through pragmatic testing and basic analysis, the intent is to gain insights of the performance of the schemes considered without insisting on the satisfaction of TVD. The reasons are twofold. No known technique is available to enforce TVD for a

23

high order time stepping scheme; furthermore, for problems with substantial source terms, the TVD scheme does not work well anyway (Hirsch 1990).

The results of the linear and nonlinear Burgers equations are presented in the following sections.

2.3.1 Linear Burgers Equation

The linear Burgers equation is chosen because it is relevant to many heat transfer applications. It is known that at steady state, both first-order and secondorder upwind solutions are wiggles-free while neither second-order central differencing nor QUICK is wiggles-free (Shyy 1985). Their relative merits in the context of the unsteady problems are assessed here. The time stepping schemes employed here are forward and backward Euler schemes. The A-B/A-M scheme will be included later in the nonlinear case.

For the linear problem, with C = 1 and P = 10, either 21 or 81 nodal points are utilized in the computations. The results with 21 nodes are first presented and compared in Figures 2-1 and 2-2 for the two time stepping schemes at t = 0.3,0.8, and 2. Using the backward Euler scheme, as shown in Figure 2-1, both the firstorder and second-order upwind solutions are wiggles-free, but the first-order upwind is a little more dissipative. For the steady state solutions of Eq. (2.1), it is well known that the central differencing scheme exhibits oscillations if the cell Peclet number, P, exceeds 2 (Fletcher 1988 and Shyy 1985). Figure 2-1 shows that this scheme also produces nonphysical oscillations behind the sharp front for the unsteady

24

equation. The QUICK scheme, which combines the characteristics of the central differencing and upwind schemes, also exhibits oscillations in the regions of sharp gradient; its solution profiles lie between those of the central differencing and upwind schemes.

Figure 2-2 compares the solutions obtained by using the backward Euler and Crank-Nicolson schemes; they are noticeably different. For the first-order upwind scheme, with C = 1, the initial profile is better preserved and the gradient is less smeared out with Crank-Nicolson. In conjunction with the first-order upwind

25

scheme, Crank-Nicolson, being less dissipative but yielding no wiggles, performs better than backward Euler. Similarly, the performance of a second-order time stepping scheme such as Crank-Nicolson is less satisfactory when combined with other second-order convection schemes due to the fact that both time stepping and convection schemes are dispersive. As illustrated in Figures 2-2 (b) -(d), nonphysical solutions are observed for all three second-order convection schemes. For both the central differencing and QUICK schemes, the oscillations yielded by Crank-Nicolson are substantially larger than those by backward Euler. For the second-order upwind scheme, while the Crank-Nicolson scheme yields improved solutions with sharper gradients and no oscillations at t = 0.8, it produces a noticeable nonphysical undershoot at t = 0.3.

The results presented in Figures 2-1 and 2-2 are generic and not affected by the variations in the number of grid points. In order to demonstrate this fact, Figure 2-3 compares the performances between the Crank-Nicolson scheme and the backward Euler scheme, along with the four different convection approximation schemes, where the values of P and C remain the same as before (10 and 1, respectively), but the number of nodal points is increased from 21 to 81. All the major characteristics depicted in Figure 2-3 can be observed in Figure 2-2.

2.3.2Nonlinear Burgers Equation

For the nonlinear Burgers equation, two different values of viscosity, P =

5x10-3 and 5x10-5, are used to explore the effect of solution gradients on the

26

performance of each scheme. Three time stepping schemes, including the first-order backward Euler scheme, the second-order Crank-Nicolson scheme, and the thirdorder Adams-Bashforth/Adams-Moulton predictor-corrector scheme, are studied along with the first-order and second-order upwind schemes. Solutions obtained at two time instants, i.e.,t = 2 and 3.5,are presented.

The case of P = 5x103 is considered first. Figure 2-4 compares the solutions yielded by the three time stepping schemes in conjunction with (a) the first-order upwind convection scheme and (b) the second-order upwind convection scheme. All the computations shown in Figure 2-4 are based on 51 nodes (Ax = 0.02), and At = 0.01, resulting in P = 4u and C = 0.5u. Figure 2-4 shows that with first-order upwind, the numerical solutions yielded by backward Euler and Crank-Nicolson are fairly comparable while noticeable improvement can be obtained by using the thirdorder A-B/A-M scheme. Figure 2-4 indicates that a first-order scheme in space and a third-order scheme in time can actually form a good combination.

Regarding the use of the second-order upwind convection scheme, shown in Figure 2-4 (b), all three time stepping schemes perform well, producing solutions in good agreement with the analytical one. In order to measure more quantitatively the performance of each individual scheme under study, Figure 2-5 depicts an L-2 error norm with respect to Ax, where the Courant number is held fixed for all cases, i.e., At and Ax vary in exact proportion. The L-2 error norm, e2(t), is defined as follows,

27

N
(0' (i't)-4(i't))2 (2.19)
e2(t) =
N

where ,4(i,t) is the exact solution, 4(i,t) is the numerical solution, and N

represents the total number of nodal points.

Several features are revealed by Figure 2-5. First, with a given convection scheme, the orders of accuracy of all three numerical solutions, with respect to Ax, yielded by the three time stepping schemes are essentially the same. While the local truncation error analysis through a standard Taylor series expansion shows that the numerical accuracy of the first-order upwind convection scheme is proportional to Ax, Figure 2-5 indicates that for the present case the errors are actually closer to

(Ax)5 instead of Ax. It is noted that the solutions shown in Figure 2-5 are based on the fixed value of C, hence both At and Ax contribute to the characteristics observed there; here the overall accuracy is higher than that indicated by a local truncation error analysis based on the spatial accuracy of the first-order upwind convection alone. As to the second-order upwind convection scheme, Figure 2-5 (b) shows that the solutions using both Crank-Nicolson and A-B/A-M yield errors proportional to

(Ax)2, while that using backward Euler is only slightly worse. Hence, for the firstorder upwind convection scheme, the resulting accuracy of all solutions are better than proportional to (Ax), while for the second-order upwind convection scheme, they are essentially proportional to (Ax)2. Overall, Figures 2-4 and 2-5 clearly

28

demonstrate that the combined performance of the discrete operators both in time and in space may not be the same as that suggested by a local truncation error analysis of the individual operator.

The case of v = 5x lO represents a problem of relatively smooth solution profiles. A more stringent case of v = 5x105 is considered next. Figure 2-6 compares the numerical solutions corresponding to the different combinations of the convection and time stepping schemes. Comparing Figures 2-4 and 2-6, it is clear that the degrees of satisfaction yielded by these schemes between the two values of v are very different. With a reduced P, the solution gradients stiffen, and the numerical wiggles are more prone to appear. Figure 2-6 (a) shows that, even with the use of the first-order upwind convection scheme, which is viewed as usually too dissipative, nonphysical wiggles develop when the A-B/A-M time stepping scheme is employed. With either backward Euler or Crank-Nicolson, on the other hand, the solutions are of smooth, diffusive characteristics that exhibit no overshoots. The difficulties of resolving the sharp gradients become more pronounced for the secondorder upwind scheme, shown in Figure 2-6 (b), in which the numerical wiggles appear in all solutions.

The L-2 error norms of the solutions of v = 5x105 are summarized in Figure 2-7 where the Courant number again remains unchanged as Ax is varied. Unlike those curves shown in Figure 2-5 for v = 5x103, the results in Figure 2-7 for v = 5x105 are of inconsistent appearances, depicting no recognizable correlations with Ax among different discretization schemes or between different time instants. In

29

particular, the solutions associated with the A-B/A-M scheme are most erratic; they can either be less accurate than the solutions using other time stepping schemes or produce more error as Ax is reduced. The results presented so far clearly suggest that the three time stepping schemes exhibit different dispersive characteristics; as the leading local truncation error of a time stepping scheme becomes higher-order in At, the resulting solutions tend to be more dispersive, in form of either long wave overshoot or short wave wiggles.

The above discussions are based on a single set of the Courant and cell Peclet numbers. Since each of the time stepping and convection schemes possesses its own diffusive and dispersive characteristics, it appears useful to adjust the values of P and C so that a time stepping scheme of particular dispersive and diffusive characteristics can be utilized in conjunction with a convection scheme of complementary characteristics to improve the overall solution accuracy. That such a possibility does exist is illustrated in Figure 2-8, where both v, 5x105, and Ax, 0.02, remain the same as those used in Figure 2-6, while At is increased from 0.01 to 0.05, resulting in a C = 2.5u. By increasing the size of At, the backward Euler scheme exhibits stronger diffusion characteristics while the other two higher-order time stepping schemes exhibit stronger dispersive characteristics. Hence, as demonstrated in Figure 2-8 the first-order upwind convection scheme, being generally more diffusive, can be used along with the Crank-Nicolson scheme, a more dispersive scheme. The net outcome is that they can produce satisfactory solutions exhibiting no spurious wiggles with a larger time step, yielding saving in the computing cost. On the other hand, the

30

second-order upwind convection scheme, being less diffusive, can apparently benefit from the backward Euler time stepping scheme, which is more diffusive. These interesting features suggest that the different characteristics in convection and time stepping schemes can indeed be selectively combined to form an overall algorithm capable of suppressing excessive numerical dispersion and/or diffusion with reasonable At.

2.4 Von Neumann and Spectral Analyses In order to gain some insight of the dissipative and dispersive characteristics of a given numerical scheme, the von Neumann stability analysis (Warming and Hyett 1974, Morton 1971) and the discrete fast Fourier transform (FFT) spectral analysis are utilized. The von Neumann stability analysis is applicable to linear, constant coefficient problems with periodic boundary conditions.

The elementary solution for linear Burgers equation is

0(t,x) = e -('' + k"ue ik.x (2.20)

where i = f-l,

km = 27r/L m is the wave number of the m-th component,

and Lm is the wavelength of the m-th component.

Based on the elementary solution, the amplification factor of the exact solution can be defined

31

C X~
G = (+_t) - e 7e -ioc (2.21)
0 (t)

where IG, = exp(-C32/p) is the amplification magnitude, 0 = - OC is the phase angle, 3 = km, Ax, and 21rAx/o is the wavelength.

It is clear that G, depends on the cell Peclet number, Courant number, and the wave number. The amplification factors yielded by the von Neumann analysis for the backward Euler and Crank-Nicolson time stepping schemes along with various convection schemes for linear Burgers equation are given in Tables I and 2. respectively.

In order to illustrate the characteristics of each scheme more clearly, a series of results are computed for P = 10', C = 0.5, Ax = 0.05, and a periodic boundary condition. The initial condition is defined as

O(x,0) = , for 0.05:5 x 0.2 o r 0.55 x 0.7 0, otherwise

Each of Figures 2-9 to 2-16 includes four parts : (i) the comparison of exact and numerical solutions at t = 2At, (ii) the spectral analysis of the numerical solutions yielded by FFT, and the comparison of (iii) the amplification factors as well as (iv) the phase angles according to the von Neumann analysis. The performance of the backward Euler time stepping scheme along with different convection schemes are presented in Figures 2-9 to 2-12; the corresponding results of the Crank-Nicolson time stepping scheme are presented in Figures 2-13 to 2-16.

32

For the linear problem subject to a periodic boundary condition, the ratios of magnitudes at each wave number between two consecutive time steps shown in the spectral plot (iii), are identical to the amplification magnitude shown in (iii).

It is noted that with a high Peclet number, P = 10', and a periodic boundary condition, the exact solution essentially maintains its initial profile just like a pure wave equation. Consequently, the spectral distributions of the exact solution remains unchanged with time,and the amplification magnitude shown in (iii) are equal to 1 at all wave numbers. Several observations can be made by inspecting Figures 2-9 to 2-16. First, Crank-Nicolson does not yield less dissipation than the backward Euler scheme for all wave numbers. In fact, except for the case of the central difference convection scheme, part (iii) of Figures 2-9 to 2-16 demonstrate that for all three upwind type of convection schemes, Crank-Nicolson is less dissipative than backward Euler only in the long wave length regime; with shorter wave, it is actually more dissipative than backward Euler. Secondly, the amplification magnitudes of the second-order central difference scheme are very close to the exact values; when combined with Crank-Nicolson, it is indistinguishable from the exact profile, i.e. there is no numerical damping created by the scheme, which is good. However, with the central difference scheme, the short waves travel with much smaller phase speeds than the exact solutions. It is these phase speed errors in the short wavelength that cause the central difference scheme to perform unsatisfactorily. This observation explains why the central difference scheme usually yields solutions with highly noticeable 2Ax oscillations.

33

The other interesting aspect revealed by the present analysis is that the second-order upwind scheme is actually more dissipative than the first-order upwind scheme at short wavelengths as indicated by part (iii) of Figures 2-9 and 2-11. With both Crank-Nicolson and backward Euler, the second-order upwind scheme is less dissipative than the first-order upwind scheme only in the long wave regime. This feature may be useful if a solution profile contains sharp gradients because the phase speed errors associated with short waves are less pronounced. Furthermore, by comparing Figures 2-11 and 2-15, it becomes clear that it is the phase errors in the long wavelength that are mostly responsible for the inaccuracies created by the combined Crank-Nicolson/second-order upwind schemes. This can explain why its solutions exhibit a single long wave overshoot instead of 2Ax short wave oscillations. For the first-order upwind scheme, it appears that in the long wave regime, CrankNicolson is not only less dissipative but also more accurate in phase speed. Accordingly, as already concluded previously, backward Euler/second-order upwind or Crank-Nicolson/first-order upwind form good combinations. The performance of QUICK lies between the central difference and upwind schemes.

2.5 Summary

Both linear and nonlinear Burgers equations have been used as the model problems to investigate the interaction of the time stepping and convection schemes for unsteady flow simulation. The von Neumann stability analysis and the FFT spectral analysis are also utilized 'to aid our understanding of the performance

34

characteristics of these schemes. Interesting features are revealed by these analysis. For example, Crank-Nicolson is less dissipative than backward Euler only in the long wave regimes; it is actually more dissipative as the wavelengths are shortened toward 2Ax. Similarly, with the same time stepping scheme, the second-order upwind scheme is less dissipative than the first-order upwind scheme mainly for the longer waves. Consequently, for the second-order upwind scheme, backward Euler yields more accurate phase speed than Crank-Nicolson; while for the first-order upwind scheme, it is Crank-Nicolson that is more satisfactory in phase speed. Regarding the second-order central difference scheme, because it contains very little numerical damping for all wavelengths, the large phase speed errors, especially at short wave lengths, cause its solutions to exhibit pronounced 2Ax short wave oscillations. In contrast, for second-order upwind, because of its higher dissipation rate at short wavelengths, it is the phase errors of longer waves that cause the noticeable inaccuracies; consequently, the solution error tends to appear as a single long wave overshoot. The performance of QUICK essentially lies between the central differencing and upwind scheme.

Overall, among all the schemes tested, either a combination of first-order upwind for convection and Crank-Nicolson for time, or a combination of secondorder upwind for convection and backward Euler for time performs better. A consistent second-order accuracy for both time and space discretizations does not produce good results for the time dependent flows with large gradients. The results suggest that a selective combination 6f the convection and time stepping schemes of appropriate dispersive and diffusive characteristics can attain good performance.

Table 2-1 Amplification Factors of Backward Euler Scheme
for Linear Burgers Equation
scheme amplification factor
first-order p
upwind scheme
G = C
[-(2+P+-) + (2+P)cos ] - i(Psin )
C

second-order p
central difference

[-(2+-) + 2cosp I - i(Psinp)
C

second-order p
upwind scheme G -[-(2+ P+-) - -cos(2p) +(2+2P)cosp I - i[-sin(2P) -2PsinPj
2 C 2 2

QUICK
scheme

[-(2+-P+-) - -cos(2) + (2+-)cosP I - i[-sin(2P) -5sinpJ
8 C 8 2 8 4

Table 2-2 Amplification Factors of Crank-Nicolson Scheme for Linear Burgers Equation

scheme amplification factor
first-order
upwind scheme [(2 +P-2-) - (2+P)cosp] + i(Psins)
G5 C
(2 +P+2-) + (2+P)cos p I - i(Psinp)
C

second-order P
central difference [(2-2-) - 2cosp] + i(Psin p)
G6 C
[-(2+-) + 2cosp I - i(Psinp)
C

second-order 3 P P P
upwind scheme [(2+P-2-) + cos(2P) -(2+2P)cosPI + i[- sin(2p)+2PsinpJ
G7 - 2 C 2 2
[-(2+ P+2 )- cos(2 P) + (2+2P)cos p] + i [sin(2p ) - 2Psinp]
2 C 2 2

QUICK 3 P P P P 5
scheme [(2+-P-2-)+- cos(2p)-(2+-)cosP] + i[- sin(2p)+ -Psinp]
[( 8 C 8 2 8 4
[(2+-3P+2) - - cos(2 P) + (2 +-P)COSj I + i[-Psin(2 P) --5sinpi
8 C 8 2 8 4

1.2

U.2 U.4 U.0 0 I.M
x

1.8

1.4 1.0 0.6

0.2

1.0 0.8 S0.6

0.4 0.2 a 0

1.2 1.0 0.8 4 0.6

0.4 0.2

2

1 3

(a) 1 = 0.3

3
12

1 first- order upwind scheme
2 second-order central differencing scheme
3 second-order upwind scheme
4 QUICK scheme

0.6 0.8
x

1.0 1.2

(c) t = 2 (steady state)

1.2

(b) t = 0.8

Figure 2-1 Comparison of solutions of unsteady linear Burgers equation using different
convection schemes, backward Euler time stepping, P = 10, C = 1, and 21
nodal points.

2

0 0.2 0.4

C

0

b

V

1.2

1.0 0.8 0.6

0.4 0.2

0

(a) first-order upwind scheme

1.2

1.0 0.8 0.6

0.4 0.2

0

0

0.2 0.4 0.6 0.8 1.0
x

(c)

a

1.2 0 0.2 0.4 0.6 0.8 1.0
x

t 03

second-order central differencing scheme d. =0.8 1=03

1.2

(b) second-order upwind scheme (d) QUICK scheme

Figure 2-2 Comparison of solutions of unsteady linear Burgers equation obtained using
the backward Euler and Crank-Nicolson time stepping schemes.

Crr

S00.3

b

1.2 1.0 0.8 0.6 0.4 0.2

0

0 0.2 0.4 0.6 0.8 1.0 1.2
x

(b) second-order upwind scheme

, =0.8

-, /t = .3
I. /

a

(a) first-order upwind scheme

0

0.2 0.4 0.6 0.8
x

1.0

1.2

(d) QUICK scheme

Figure 2-3 Comparison of the backward Euler and Crank-Nicolson time stepping
with four convection schemes for unsteady linear Burgers equation
nodal points, and P = 10, C = 1.

schemes using 81

1-08

izO
c ~-

(c) second-order central differencing scheme

1.2 1.0 0.8 0.6 0.4 0.2

0

t=0.3

o0

t=08

d

0.4 0.3 D 0.2

0.1

0

0.4 0.3

p 0.2

0.1

0

0 0.2 0.4 0.6 0.8
x

1.0 1.2 0
b

Lxact t 2
solution
3

2

= 3.5 Exact
solution 3

2

L .---- --J

0.2 0.4 0.6
x

(a) first-order upwind convection scheme

(b) second-order upwind convection scheme

Figure 2-4 Effect of time stepping schemes on solution accuracy of the nonlinear
Burgers equation with v = 5x103 , Ax = 0.02,and At = 0.01.

1 2

2
3

1 Backward Euler
2 Cronk -Nicolson t 3 5
3 A-B/A-M

I i iI

a

0

0.8 1.0 1.2

10.2

I 4

10-i 102 10-

t 2

2 3

1 Backward Euler 2 Crank -Nicolson
3 A-B/A-M

= 3.5

2

,ag3

104 [

10-5
10-1 10
Ax

t=2

2
3

t =3.s

3

1 102 10
b Ax

(a) first-order upwind convection scheme

(b) second-order upwind convection scheme

Figure 2-5 The L-2 error norm with respect to Ax of the nonlinear Burgers equation with
= 5x103 and At/Ax = 0.5.

a

-1

0

0 0.2 0.4 0.6
X

0.8 1.0 1.2 0
b

1 Backward Euler 3 t=2
2 Crank -Nicolson
3 A-8/A-M Exact
solution
1 2

0.2 0.4 0.6 0.8 1.0
X

(a) first-order upwind convection scheme

(b) second-order upwind convection scheme

Figure 2-6 Effect of time stepping schemes on solution accuracy of the nonlinear
Burgers equation with v = 5x105 , Ax = 0.02, and At = 0.01.

0.4 0.3

p 0.2

0.1

0

0.4 0.3

lz2

3
2

It 3.5

2

0

a

t 3.5

3
Nxcl
2 ,slto

t'~)

1.2

-_ 0.2[

0.1[

102 10

10. 1t.-

Ax

Cj.)

101

Ax

10-1 10.2
b

(a) first-order upwind convection scheme

Figure 2-7

(b) second-order upwind convection scheme

The L-2 error norm with respect to Ax of the nonlinear Burgers equation with v = 5x105 and At/Ax = 0.5.

- 3t=2
-3
2

eo

t 2

3

1 Backward Euler 2 Crank-Nicolson
3 A-B/A-M

3 t =3.5

2

1010 10.2

a

3~ 1 3 5
3~3
2

0.4 0.3 D 0.2

0.11

0

0 0.2 0.4 0.6 0.8 1.0 .
X

0.40

0.36 P 0.32

0.28

0.24

a

0.60

X

0.65

0.75

,-- - - - t2

0.60

b

X

0.65

0.75

(a) first-order upwind convection scheme (b) second-order upwind convection scheme

Figure 2-8 Effect of time stepping schemes on solution accuracy of the nonlinear
Burgers equation with P = 5x10- , Ax = 0.02, and at = 0.05.

I

1 Backward Euler 2 Crank - Nicolson

t= 2

3A-B/A-M

Exact solution

. -. t 2

0 0.2 0.4 0.6 0.8 1.0 1.
X

1=2

Exact solution

3 /
- 2

1.5

1.0

- 0.5

0

P = o' C = 0.5 Ax 0.05

=2Ai

I3

-0.5 - I
0 0.2 0.4 06 0.8 1.0 1.2
x
2x/P(=wavelength/Ax)

(a) comparison of exact solution and
numerical solution at t = 2 At

1.2 1.0 0.8 p0.6

0.4 0.2

0

0

1 2 3 4 5 6 7 8 9 10
2x/p(= wavelenigIh/Ax)

(c) amplification magnitude

0.50
0.45 (a)S=0
0.40 (b) I = At
(c) S = 2AI
0.35
0.30 - exact
0.25 - - - numerical
0.20
0.15
0.10 (b)
00

0 1 2 3 4 5 6 7 8 9 i

)

2n/P(=waveength/Ax)

(b) spectral analysis at several
time instants

0

-0.5

-10

-1.5

-2.0
0

(a)

(a) C = 0.25
(b) C - 0.5

1 2 3 4 5 6 7 8 9 10
2i/P( =wavele ngIi/Ax)

(d) phase angle

Figure 2-9 Performance characteristics of the combined backward Euler and first-order
upwind schemes with periodic boundary condition, P = 10' , C = 0.5,and
Ax = 0.05.

, .-- -- - - - - -

(b)

(.) C = 0.25 (b) C = 0.S

LA

1.5 1.0

- 0.5

0

-0.5
-a

0.8

1.0 1.

4A

2

2n/p(=wavelength/Ax)

(a) comparison of exact solution and
numerical solution at t = 2 At

1 2 3 4 5 6 7 8 9
2K/P(=wavelenglh/Ax)

(c) amplification magnitude

0.5

0

-0.5

- 1.0

- 1.5 t

- 2.0 t 0 0

(a) I = 0
(b) I = at
(c) I = 2At

P t0'
-5 C =0.5
Ax = 0.05
F .1 1=2AI
-5 5 .5

* . - .',,

cal

123. 5678

1 2 3 4 5 6 7 8 9
2x/p(=wavelengIh/Ax)

(d) phase angle

I
10

Figure 2-10

Performance characteristics central difference schemes C = 0.5,and Ax = 0.05.

of the combined backward Euler and second-order with periodic boundary condition, P1 = 10' ,

0.2 0.4 0.6
x

u.ZU
0.45
0.40 0.35 0.30 0.25
0.20 0.15 0.10 0.05
0 0

1.2 1.0 0.8 1 0.6

0.4 0.2

0

Ill

---------- (b).-----(a) C - 0.25
(b) C = 0.5

U

Iv

(b

(a) C - 0.25
(b) C - 0.5

C\

1 2 3 4 5 6 7 8 9 10

2n/P(=wavelength/Ax)
(b) spectral analysis at several
time instants

- exact
- - - nunerli

(b)

aI

I

1.5

1.

0.

0.

0

0

0

5

0

1.

*0.

-0.

1.

Figure 2-11

P = to, C = 0.5 Ax = 0.05 t=2AI

0

-0.5

- 1.0

-1.5

' 0.2 0.4 06 0.8 1.0 1.2 2x/P(=waelengh/Ax)

(a) comparison of exact solution and numerical solution at t = 2 At

2
Ill

0 ,-- ' . --- -

6 ,.(
- (b)
.4,-'

.2 (a) C = 0.25
(b) C = 0.5

0 1 2 3 4 5 6 7 8 9 10 2n/P(=wavelenglh/A x)

(c) amplification magnitude

1 2 3 4 5 6 7 8 9
2n/p(=wavelength/Ax)

(d) phase angle

-4

10

Performance characteristics of the combined backward Euler and second-order upwind schemes with periodic boundary condition, P = 10' , C = 0.5, and Ax = 0.05.

- 2.0 L
0

Iv

(a) (b) ---(b)

(a) C = 0.25
(b) C 0.5

0.50
0.45 0(a)=0
(b) I = At
0.40 (c) t =2At
0.35
0.30 - exact
0.25 .--- numerical
0.20
0.15
0.10 (b)
0.05 (C

00 1 2 3 4 5 6 7 8 9 10

2n/p(=wavelength/Ax)
(b) spectral analysis at several
time instants
0.5

0.2 0.4 06
x

P = o,' C =0.5
-'I~A /' 0.05
I= 2At

0.8 10 1.2

2n/P(=wavelength/Ax)

(a) comparison of exact solution and numerical solution at t = 2 At

2
i
0 -. ----------- :::::::::::
(a) - ---(b)
.6

4

.2 (a) C - 0.25
(b) C a 0.5
0
0 1 2 3 4 5 6 7 8 9 1 2n/P(=wavelenvglh/Ax)

(c) amplification magnitude

1.5 10

IV

0

-0.5

- 1.0

0.50
0.45 0.40 0.35 0.30 0.25 0.20 0.15 0.10 0.05
0

1 2 3 4 5 6 7 8 9
2n/P(=wavelength/Ax)

(d) phase angle

00

10

Figure 2-12 Performance characteristics of the combined backward Euler and second-order
QUICK schemes with periodic boundary condition, P = 103 , C = 0.5, and
Ax = 0.05.

4,
-a

a.

1 I 2 3 4 5 6 7 8 9 10

2n/P(=wavelength/Ax)

(b) spectral analysis at several
time instants

e 0.5

0

-0.5
0

II (a) I = 0

(a) I = 0
(b) = A I
(c) t=2AI

- exact
- - - nunerli

(a)
(b) 'e.

1 .

1.

0.

0

0.

(b)

(a) C = 0.2S (b) C - 0.5

-

a

Cal

05

-1.0

0 0

1.5 1.0 +0,5

0

-

0 0.2 0.4 0.6 0.8 1.0 1.2
x
2n/O(= wavelength/Ax)

(a) comparison of exact solution and
numerical solution at t = 2 At

1.2 1.0 0.8 0.6

0.4 02

(@) C =0.25
0 (b) C = 0.5
0 1 2 3 4 5 6 7 8 9 10 2I/p(=wavelengih/Ax)

(c) amplification magnitude

0.50
0.45 0.40 0.35 0.30 0.25
0.20 0.15 0.10 0.05

0

0.5

0

-0.5

- 1.0

-1.5

- 2.0

P = l0' C = 0.5 . = 0.05 I 2AI

0 1 2 3 4 5 6 7 8 9 10
2n/P(=waveleigIh/Ax)

(d) phase angle

Figure 2-13 Performance characteristics of the combined Crank-Nicolson and first-order
upwind schemes with periodic boundary condition, P = 10' , C = 0.5,and
Ax = 0.05.

- exact
- - - numerical

(b)

1 2 3 4 5 6 7 8 9 1

2n/P(=waveengIh/Ax)

(b) spectral analysis at several
time instants

-(b)

(b)

(a) C - 0.25
(b) C = 0.5

II (.) 1=0

(a) I = 0
(b) I =At
(C) I 2AI

-
.- -- ---------
(b)

4V

HI

0.50
0.45 0.40 0.35 0.30 0.25
0.20 0.15 0.10 0.05

0.2 0.4 0.6 0.8 1.0 1.2 U0 1 2 3 4 5 6 7 8 9 10

2n/P(=waveIengIh/Ax)

(a) comparison of exact solution and
numerical solution at I = 2 At

0 1 2 3 4 5 6 7 8 9
2 n/P( =wavelengIh/Ax)

(c) amplification magnitude

2n/P(=wavelength/Ax)

(b) spectral analysis at several
time instants

CD

10

0.5

0

-0.5(a - ---(b---- - -

- 1.0 Mb

-1.5 (a) C - 0.25
(b) C - 0.5
- 2.0
0 1 2 3 4 5 6 7 8 9 10
2n/P(=wavelength/Ax)

(d) phase angle

Figure 2-14

Performance characteristics of the combined Crank--Nicolson and secondorder central difference schemes with periodic boundary condition, 11 = 10', C = 0.5, and Ax = 0.05.

1.5

1.01

-- 0.5

0

-0.5

.( .' P = to,'
!! II C =0.5
Ax = 0.05
.- .- I = 2

* .

(a) = 0 (b) I = At (c) I =2At

- exact
- - - numerical

(b)

1 2

1.0 0.8 p0.6

0.2

0

(a) (b)

(a) C = 0.25 (b) C - 0.5

LJ~
0

0

1.5 1.0 + 0.5

0

x
2n/P(=waveIeng(h/Ax)

(a) comparison of exact solution and
numerical solution at t = 2 At

1 2 3 4 5 6 7 8 9
2n/P(=wavelengIh/ax)

(c) amplification magnitude

P - to' C = 0.5
-. fI= 0.05
I . 2A 1

0 0.2 0.4 0.6 0.8 1.0 1.

2

0.50
0.45 0.40 0.35 0.30 0.25
0.20 0.15 0.10 0.05
0

1 2 3 4 5 6 7 8 9

2n/P(=wavelength/Ax)

(b) spectral analysis at several
time instants

0.5

0

0.5(b . 1.0 b)

-1.5 (a) C 0 .25
(b) C 0.5

10

0

1 2 3 4 5 6 7 8 9
2n/P(=wavelength/Ax)

)

10

(d) phase angle

Figure 2-15

Performance characteristics order upwind schemes with and Ax = 0.05.

of the combined Crank-Nicolson and secondperiodic boundary condition, P = 10' , C = 0.5,

(b)
(b)
(c)

(b) t = At
(c) I = 2Ag

- exact
- numerical

4J

S
N

-0.5

1.2 1.0 0.8 S0.6

0.4 0.2

n

Il

(b)

(a) C - 0.25
(b) C = 0.5

0

-

1.5,

1.01 + 0.5

0

-0.5l
0 0.2 0.4 0.6 [.s 1.0 1.

1.2 1.0

0.8 p0. 6

0.4 0.2

0.50
0.45 0.40
0.35 0.30 0.25
0.20
0.15 0.10 0.05

2

x
2n/P(=wavelenglh/Ax)
(a) comparison of exact solution and
numerical solution at t = 2 At

0 1 2 3 4 5 6 7 8 9
2X/p(=wavelength/Aax)

(c) amplification magnitude

U0 I 2 3 4 5 6 7 8 9 10

2n/P(=wavelength/Ax)
(b) spectral analysis at several
time instants

0

-0.5

- 1.0

-1.5

-2.0
10 0

(b)

(a) C - 0.25
(b) C a 0.5

1 2 3 4 5 6 7 8 9 10
2n/P(=wavelength/ Ax)

(d) phase angle

Figure 2-16 Performance characteristics of the combined Crank-Nicolson and second-order
QUICK schemes with periodic boundary condition, P = 10' , C = 0.5, and
Ax = 0.05.

10' 0.5 = 0.05
2A(

II (a I = -

(a) 1 0
(b) I At
(c) i 2AA

- exact
- - - numerical

(a) (b,
A\ (C)

(a."

(b)

(a) C - 0.25 (b) C = 0.5

- . . . . . . .

I

U

P .

C .
, All

CHAPTER 3
MENISCUS CHARACTERISTICS WITH APPLICATION TO THE EFG PROCESS

3. 1 Introduction

As already described, during the course of fibre growth, a meniscus is formed adjacent to the solidified material. In order to obtain crystals of definite diameters with given growth rates or rates of pulling, the control of the behavior of the meniscus of liquid bridging the elevated solid and the melt plays a key role. In this chapter, we intend to examine the behavior of the meniscus and shed some light on fundamental aspects of generic interest that have hitherto received only diffuse attention.

Studies of meniscus profiles and properties have been made by various investigators. In the works of Katoh et al. (1990), and Pitts (1976), and Huh and Scriven (1960), the meniscus is formed by dipping solids of different geometries into baths of infinite extent. On the other hand, in the crystal-pulling literature, theoretical treatments of meniscus shape and stability abound (Michel 1981, Surek et al. 1980, Tatarchenko and Brenner 1980), mainly in the form of linear analysis, and simulations of crystal growth processes including heat transfer have been made (Backman 1992, Thomas et al. 1986), with the consideration of equilibrium meniscus shape. However, the physics involved in this process, in particular, the

53

54

thermodynamics governing the menisci and the issues related to contact angles have not been adequately clarified (Dussan V. 1979). Furthermore, to our knowledge, the possibility of multiple meniscus shapes for a single contact angle and consequences thereof has not been studied. It is generally accepted (although not without argument) that both from considerations of force and energy, the equilibrium contact angle may be obtained from the well known Young's formula.

Under debate is the concept of the dynamic contact angle, and that of advancing and receding contact angles (Miller 1985, Ngan and Dussan V. 1982, Schwartz and Tejada 1982). Dynamic contact angles are a feature of concern since in such processes the solid to which the meniscus is connected is subjected to pulling as well as rotation while being retracted from the melt. Specifically, one seeks to determine a functional form for the variation of contact angle 0, with the velocity of the contact line. While success has eluded theoretical treatment, experiments have provided a broad picture of the behavior of the apparent contact angles with contact line motion. Hysteresis of dynamic contact angles implies that the behavior of Oc depends on the direction of motion of the contact line. The static contact angle achieved at a contact line that resulted from advancing fluid, say 0, is greater than that formed by receding liquid, say OR. Thus, the apparent angle assumed under static condition by the liquid-gas surface at the trijunction is either 0, or OR depending on how it was established. When the contact line is compelled to move, it will remain pinned while the contact angle ( via a series of jumps) traverses the gap 01 to OR or vice-versa, depending on the direction of impending motion.

55

Thereupon, the contact line will move, and a slip condition (de Gennes 1991, Koplik et al. 1989) and caterpillar motion are believed to be the characteristics that accompany the translation of the contact line. It is generally believed that the apparent contact angle monotonically increases above A for advancing contact lines and decreases below 4>R for receding contact lines, although the physics of the motion or an appropriate model (Joanny and de Gennes 1984, Jansons 1986) is not established.

In the following we will draw attention to some of the above issues concerning the meniscus formation caused by a pulling process. Although satisfactory solutions to these theoretical weak spots are awaited, an attempt is made to cast them in perspective. Herein, we once again solve the Laplace-Young equation for the equilibrium meniscus profiles. The boundary conditions are chosen to correspond to that of the edge-defined film-fed growth (EFG) process. We seek to identify the forms and properties of menisci and the parameters affecting shape and stability. In particular, the Bond number, contact condition, aspect ratio (meniscus height/base radius) and the degree of pressurization (interacting with gravity) significantly influence meniscus behavior. Our motive is to understand the equilibrium behavior of menisci from a fundamental viewpoint. As will be demonstrated, the issue of equilibrium meniscus behavior, even without consideration of dynamic contact angles, holds many points in need of elucidation.

56

3.2 Configuration, Assumptions and Formulation

The configuration chosen here is shown in Figure 3-1, where h,, r,, and rb are the meniscus height, radius of wetted region at top surface, and die radius respectively. We are particularly interested in the axisymmetric situation related to forming thin fibers (a few mm or smaller in diameter). As will be discussed, the diameter of the fibre can be quantified in relation to a characteristic capillary length scale of a given material through the Bond number. Since the geometry is

axisymmetric, the appropriate Laplace-Young equation for the liquid-gas interface is applied, which in dimensional form is (Surek et al. 1980, Pitts 1976, and Huh and Scriven 1960)

Yig 1
ylg fz 3
(1+"Of ) f (1+ )2 (3.1)

= Apgz - Ap

where f is the radial co-ordinate defining the surface of the meniscus, yi is the surface tension of the liquid-gas interface, Ap is the difference in density between liquid and gas phases, and Ap is the difference in applied pressure between the liquid, p1 , and gas phases, pg , i.e., Ap = p, - pg . In our results, below, we will express Ap as a percentage of ambient pressure, considered fixed at 1 atmosphere. Consider a capillary length defined as

57

a= __ (3.2)
\A pg

where, for example, for sapphire (A1203), yg = 0.69 N/m, Ap = 3.4x10' kg/M3 , and hence the capillary length is about 4x10-3 m. The Bond number, which relates the radius of the base or die to the capillary length, and controls the relative importance of gravity to capillary forces is given by

2
Bo=- (3.3)
a2

We specify the aspect ratio, h/rb and Bond number, Bo, in each case. The meniscus height is calculated from the given aspect ratio and re. We non-dimensionalize all physical quantities as follows, *'s representing non-dimensional values, z*= z/a , f= f/a , Ap*= Ap/(Apga), Ap * = 1. Again, using the sapphire as an example, the reference pressure p, = Ap g a is about 130 N/m2 Thus, the non-dimensionalized profile equation, after having dropped all *'s is,

____ 1
3 Z-Ap (3.4)

(1+f2 )2 f (1+f2)2

This equation is seemingly invariant with respect to Bond number, since the Bond number does not explicitly appear in Eq. (3.4). This is a misleading impression, since the Bond number does influence meniscus profiles via the term z, which is the non-dimensionalized gravity term. For a given aspect ratio, h/rb, the effect of gravity

58

increases with Bo, since a higher Bond number implies a taller meniscus. When Ap 0, the effect of Bond number will influence menisci through the relative values of hydrostatic pressure, Apgz, and externally imposed pressure difference, Ap, as will be shown from our calculations.

The Young's contact condition at the top of the meniscus is,

Y9coso = YS - Y1, (3.5)

where -y,g, -yl and yig are respectively the surface tensions of the solid-gas, solid-liquid and liquid-gas interface. , is the static contact angle measured with respect to the liquid-solid interface.

It may be noted that this relation is strictly applicable to a perfectly smooth surface without energy heterogeneity. It merely defines the static contact condition at the three-phase line for an ideal surface. However, if the solid surface is taken to consist of 'patches' of ideal smooth surfaces, we may impose local Young's contact angle at the contact line. This situation has been shown previously to offer one plausible model for hysteresis (Huh and Mason 1977, Eick et al. 1975, Neumann and Good 1972). It is also evident that such heterogeneities or those associated with crystalline orientations or defects, may cause the meniscus to deviate from axisymmetry. Contortions of the contact line are developed in order to ensure energy minimization, by locally enforcing on the liquid-gas surface the Young's angle at the contact line. This in turn distorts the meniscus surface which, for equilibrium, has to conform locally to the Young's condition. In our work, we avoid the complications associated with these aspects by assuming an ideal surface at the top. It is a

59

conventional practice to account for the contact angle hysteresis by adopting a relation such as Wenzel's, namely

cos4, = r, cos4 (3.6)

where rf is a roughness factor and O.,, is the apparent contact angle (Oliver et al. 1977). So as not to dilute our focus on the fundamental meniscus behavior under equilibrium conditions, this practice is not adopted. It is possible, however, by scanning through the available range of 4 for solutions to the Laplace-Young equation. to incorporate the effect of Eq. (3.6) by resorting to the appropriate value of 0,pp. Hence, a dynamic process with varying 4.P, can be treated simply as a juxtaposition of equilibrium solutions with varying 4,.

The Laplace-Young equation being an ordinary differential equation of second order, two boundary conditions are necessary. For the EFG process, the radius is fixed at the base. The liquid-vapor line is assumed to be pinned at the edge of the die as shown (Surek et al. 1980, Oliver et al. 1977), such pinning having been observed in practice.

In the present simplified framework, i.e.,without explicitly accounting for heat flow and solidification, we choose to picture the meniscus as attached to a flat plate. Thus, at the top, for a given meniscus height, h, , the contact angle is specified and the radius of the trijunction point can be computed. This is in cognizance of the fact that the Young's contact angle is a material property, independent of the geometry of the meniscus. Such a boundary treatment also offers a simple way to simulate fibre diameter variations. However, the fact that a Neumann boundary condition is

60

specified at the top results in the existence of multiple solutions. In fact, it is known that a boundary value problem generally does not guarantee a mathematically unique solution. Hence, guidance will be needed to help select the solution that is physically stable. The implications of mathematically non-unique solutions on the realizability of the meniscus shape will be discussed later.

3.3 Energy Formulation

The thermodynamics of the meniscus may be studied by defining the energy change in forming the meniscus. The total energy E of this system is given by the sum of the gravitational potential energy and the free energy of the surface. To nondimensionalize energy,

E* E (3.7)
A pga4

Total energy resulting from the formation of the meniscus, E, is equal to the potential energy (from the effective head) E, plus surface energy of meniscus E. plus surface energy of the wetted area of the top plate E3. The first component of the energy, that due to gravity and imposed pressure is, El = fPf2 (Apgz+Ap) dZ (3.8)

The second component, corresponding to that expended in forming the free surface is,

61

E2 = f2n yIgf(1+fz)11 dz (3.9)

The third component, the energy of wetting of top surface is, E3 = (2y-ys) r (3.10)

where r, is the radial location of the trijunction at the top of the meniscus.

We follow the usual practice of replacing the difference, (-yl, -'s) in the third component using the Young's contact condition (Pitts 1976, and Huh and Scriven 1960). This procedure facilitates evaluation of this component, since the values of contact angle 00 and -yg are often the only experimentally obtainable quantities. Thus, the third component of energy may be written as

E3 = Yig cos*, n r, (3.11)

The total energy then is, after nondimensionalizing and dropping all *'s, E = f.f(z+Ap)dz
(3.12)
+ f2 nft1 +f2)1/2dz + Irrc cos6,

By scanning through the range of possible solutions the total energy E can be computed as a function of 0,, the contact angle which is contained in the solution for meniscus shape, including that corresponding to the specified value 0.. The decision regarding stability of the meniscus can then be arrived at based on whether the specified value of 'k0 minimizes E or not.

62

3.4 Numerics

The Differential equation governing the meniscus shape, the Laplace-Young equation, is nonlinear and has to be solved numerically. The integration of the nonlinear ordinary differential Eq. (3.4) is performed using the fourth-order Adams predictor-corrector method, using the Runge-Kutta method for starting. Onehundred and one points were used along the z- direction, as deemed sufficient from grid refinement studies. The procedure adopted was to use the shooting method beginning from the fixed base at the edge of the die. Shooting was performed at various starting angles, henceforth

3.5 Results and Discussion

As previously mentioned, the meniscus profiles are significantly influenced by Bond number (Bo), applied pressure (Ap), aspect ratio (hr /rb), and the value of contact angle. Our procedure was to specify the Bond number and aspect ratio, thereby obtaining the value of rb from Eq. (3.3) , and meniscus height from aspect ratio and rb.

As already mentioned, an interesting feature of our solutions for this model of Edge-defined Film-fed growth is the presence of multiple solutions for a given

63

contact angle. Unless a methodical sorting is performed, a numerical solution procedure may not pick up all the possible solutions. In the work of Katoh et al. (1990), multiple solutions seem to have been obtained only in certain cases. This may be due to their operating regime of parameters, such as Bond number, aspect ratio etc. From our solutions it appears that the tendency to obtain multiple solutions and interesting profile shapes is greater for taller menisci, and is also influenced by the applied pressure.

3.5.1 Basic mMeniscus Behavior with Upward Pulling

Consider a typical case (Case 1) shown in Figure 3-2. The parameters are Bo = 4x10 -, he/rb (aspect ratio) = 0.5, Ap = 0. The curves corresponding to 4 against b,, and r/rb against 4, clearly show the multiple solution behavior. Thus, there is only an interval of possible contact angles 4, that can be achieved at the top of the meniscus. Outside this range of angles no solutions exist. This is not due to instability of possible equilibrium profiles, but because it is not even possible to obtain profiles satisfying the Laplace-Young equation. Also shown in Figure 3-2 are sample meniscus profiles, corresponding to 0(, of 60* and 900, respectively.

As shown in Figure 3-2, when the value of 0. is specified, and E is calculated, distinct extrema are obtained at specified values of 0.. Among the multiple solutions, the one corresponding to a minimum on the curve is the physically realizable stable meniscus profile. When either a maximum or a non-extremum point arises, we may surmise that such solutions belong to the unstable branch.

64

That the Young's condition is in fact the condition for stability is straightforward to show, say in the case of a sessile drop. In such a (constant volume) case, when gravity is ignored, the drop assumes the shape of the arc of a circle (Dyson 1978). Then the minimization result becomes merely a geometric condition compatible with minimum surface area. When gravity is considered the drop no longer remains the element of an arc , since minimization of surface free energy is not sufficient, it being necessary to minimize the total energy including gravitational potential. In the case of menisci above infinite baths, Pitts (1976) and others have shown the minimization principle via variational methods. It is the conventionally held view that Young's contact condition can be obtained both from consideration of forces as well as energy (however, see Johnson (1959) and discussion therein for criticism of this view). However, it is interesting to see that in the present configuration, multiple shapes corresponding to energy minimum, maximum and even non-extremum may all be observed at the given contact angle 0.. Furthermore, solutions to the Laplace-Young equation, which is after all an equilibrium relation across the interface, are all equilibrium solutions. From the standpoint of local energy minimization, there may be more than one stable shape. Thus the end condition and the path of development of the meniscus selects the stable solution. However, in a process in which the aspect ratio is varied, an unstable meniscus may form first. This can happen if the new trijunction location corresponding to this shape is closer to the initial trijunction location. Then, in achieving the stable solution profile for the given height, one may pass through an unstable solution profile.

65

As discussed in connection with the non-dimensional form of the LaplaceYoung equation, for low Bond numbers, the term corresponding to gravity is small compared to the curvature terms on the left hand side of Eq. (3.1). Thus the Bond number serves principally to relate the dimensions of the meniscus to the capillary length scale according to Eq. (3.3). However, for higher Bond numbers, the hydrostatic term on the right hand side of Eq. (3.1) becomes comparable to the curvature terms on the left hand side, and gravity does significantly impress on meniscus profiles, as will be illustrated in the results to follow.

For taller menisci, i.e., higher aspect ratios, solutions to menisci appear to offer a wider variety of shapes and range of angles as well as a greater tendency to multiple solutions. Retaining Bo=4x105 and Ap=O, as in Case 1, we can increase the height of the meniscus. For example, for h,/rb =2, which represents a taller meniscus, the available range of contact angles will shrink. Thus, as expected, when the height of the meniscus is increased, in the absence of pressurization, a large number of solutions to the Laplace-Young equation fail to exist.

Cases 2 and 3 (Figs. 3-3 and 3-4 respectively) are illustrative of the effects of pressure and Bond number on menisci with a higher aspect ratio where solutions with selected contact angles are shown. For Case 2, Bo=2.5x10-6,Ap=5%, he/r. =2. For this low Bond number case, from Figure 3-3, we notice that for a contact angle of 160, three equilibrium profiles are obtained. The possible contact angles at the top were in the range of 180 to 260. These profiles may not be seen in a physical situation, since the contact angles corresponding to the materials of the system may

66

not lie in this range at all. Thus, ultimately the existence of menisci will be dictated by the range of contact angles 4c at the top permitted by the materials, the range of Ob for which the pinning condition is maintained at the edge of the die, and the stability of the equilibrium solution (of the Laplace-Young equation) obtained.

The contrasting situation to Case 2 emerges in Case 3 (Figure 3-4). Here the Bo = 2.5x103(h./rb = 2, and Ap = 5% as in Case 2) is 103 times that of Case 2. The available range of contact angles is considerably widened, and a greater variety of meniscus profiles than for the lower Bo in Case 2 can be obtained. Comparison of profiles with those in Figure 3-3 shows the difference that Bond number makes to the effect of pressure. Changes in curvature arise in the profile. This implies that for the given parameters, on the right hand side of Eq. (3. 1), terms of pressure and gravity must compete to determine the shape of the meniscus. Thus, as mentioned before, for Ap 0, the Bond number can significantly influence meniscus shapes and permissible solutions. This can be seen from the scaling that we have employed in the non-dimensionalization. With fixed h, /rb, Eq. (3.3) shows that the height of the meniscus is scaled with respect to Bond number as h. oc 4Bo. Hence, for given Ap, g and a, the meniscus height increases with Bond number. For higher Bo, by definition, the relative importance of gravity (Apgy) becomes more significant compared with the pressure term. Thus, the competition between gravity and pressure in the force balance leads to a variety of curvatures for higher Bond numbers. For lower Bo, a relatively small Ap can overwhelm the gravity term, and the curvature is less sensitive to location along the z-axis.

67

Another aspect to note in Figs. 3-3 and 3-4 is that in the E vs. 4, curves both maxima and minima are obtained at specified 0. As a contrast, however, consider Case 4 (Figure 3-5). Here, we have Bo= 4x104, Ap=5% and with the hc /rb increased to 5. While the other properties appear to follow the general trend of Case 3, the behavior of energy is distinctive. Here, for 0. = 900, both minimum and maximum are obtained. However, the minimum energy state so obtained is only local. There appears to be a global minimum at 130*. For 46, = 1300, Figures 3-5

(c) and (d), the local minimum occurs at around 1000 and there appears to be a nonextremum point corresponding to the specified 4) indicating that although only one meniscus profile exists, it is not a stable one. Thus, for this case the decision regarding the stable meniscus profile becomes more involved.

In summary, the numerical simulation reveals that depending on the range of parameters a variety of solutions are possible for a given contact angle. Multiple energy minima, implying multiple stable shapes, or no minima, implying unstable shapes may be obtained for a given contact angle.

3.5.2 Effect of Pull Direction on Meniscus Behavior

The issues of the selection of a stable meniscus shape out of the possible solutions, and the implications of variation of h/rb, as well as the direction of pulling on meniscus stability, can be better addressed through the energy considerations stated earlier. In studying the variation of meniscus shape and top radius with height the effect of the following factors is investigated. They are upward/downward pulling

68

with respect to direction of gravity, Bond number (Bo), pressure difference between the liquid and gas phases (Ap), and aspect ratio (he/rb).

The results to be discussed in the following are obtained by fixing rb and varying hC/rb from 0.5 to 2.5 upward and downward, respectively. Solutions to the Laplace-Young equation are obtained for each value of h/rb and the range of available solutions in terms of 0, and resulting radius at the top (expressed as rc/rb) are determined. Figure 3-6 compares the characteristics of the meniscus for the upward and downward pulling situations, with Bo = 10' and Ap = 0. The figures display the solution range of 0, versus hc/rb, the solution range of r/rb versus h,/rb, and possible meniscus profiles for a prescribed contact angle , = 100. As shown, in the range of 4 for which the Laplace-Young equation does offer solutions, the possible maximum contact angle 0, decreases with increasing hc/rb, while the possible minimum contact angle 0c changes little during the pulling. Consequently, the solution range of 0, shrinks when h/rb is increased, implying that unless the materials comprising the system are such that the contact angle at the top has a low value, the meniscus detaches from the solid surface when he/rb is increased.

It is interesting to note that with no pressurization, while the direction of pulling causes little difference in permissible range of contact angles, it affects the trijunction location in a more profound manner. The range of permissible trijunction locations is much narrower for upward pulling than for downward pulling. This feature is illustrated by the selected meniscus profiles corresponding to 0c = 100, shown in the bottom of Figure 3-6.

69

With Bo = 10-2, when hC/rb is less than 1, the hydrostatic pressure is small relative to the curvature terms on the right hand side of the Laplace-Young equation. Hence the meniscus profiles are essentially determined from the geometrical constraints alone, meaning that within this regime the specification of r, at the bottom and 0, at the top dictates the meniscus shape. Figure 3-6 clearly indicates that the effect of pulling direction is noticeable only when the aspect ratio is high enough. The critical value of he/rb for the hydrostatic pressure to exert a significant effect is Bond number dependent. It appears that the more appropriate indicator of the importance of hydrostatic pressure is (h/rb)2Bo, i.e.,an effective Bond number defined according to h, instead of rb.

We next focus on the interaction of the direction of pulling, level of pressurization and variation of he/rb. Figure 3-7 depicts the solution characteristics for the case of Bo = 10-2, Ap = I%, i.e., there is an excess pressure inside the meniscus. Regarding the overall range of permissible meniscus solutions, both in terms of 0c and re, the direction of pulling does not make any noticeable difference, seemingly indicating that the externally imposed pressurization dominates the hydrostatic pressure. However, depending on the value of hc/rb the detailed meniscus shapes can be significantly affected by the pulling direction. For example, with hC/rb less than or equal to 1, little difference is found between meniscus shapes formed by upward and downward pulling. For h,/rb of 2 or higher, however, the differences can be substantial. Specifically, for hC/rb = 2, with upward pulling, three possible meniscus shapes exist, while with downward pulling, only one profile can be obtained.

70

When he/rb is raised to 2.5, for either direction of pulling only one meniscus profile is obtained, but the location of trijunction point for the downward pulling case differs from that of the upward pulling case.

The energy profiles corresponding to the cases in Figure 3-7 are presented in Figure 3-8; for all the values of h, /rb considered here, both directions of pulling yield energy profiles that appear essentially identical. The energy curve for h, /rb = 2 shows that for the upward pulling cases, all three possible meniscus shapes contain virtually the same amount of energy despite their noticeable differences in trijunction location. Significantly, with hc/rb of 2 or higher, the energy content for both upward and downward pulling cases at 40=70* correspond to either local maxima or nonextrema. Thus, for either direction of pulling, the meniscus shape can no longer be stable for such high aspect ratios. For larger he/rb, the E value at the specified 0 is a local maximum of the energy profile. Hence, all such menisci are not stable to small disturbances, and tend to depart from the given meniscus shape and contact angle. Since here the static contact angle does not provide a minimum energy state, the meniscus, assuming it does not altogether collapse, will tend to seek such a state by changing shape and contact angle in a dynamic manner.

3.5.3Ouasi-equilibrium Motion with Hysteresis

Contact angle hysteresis is generally acknowledged to be a consequence of three factors: (1) surface inhomogeneity, (2) surface roughness, and (3) impurities on the surface (Miller 1985). In this third set of simulations we investigate the quasi-

71

equilibrium path negotiated by the meniscus. by incorporating a simplified hysteresis model at the top trijunction, i.e., on our modeled flat plate. Specifically, the top plate is oscillated about a mean height h,' and the assumption is made that the time scale of these oscillations is such that a succession of equilibrium states is achieved. In other words, the Laplace-Young equation is solved for each position hc(t) of the oscillating plate. The contact line at the die is allowed to remain pinned at the edge, as long as the Gibbs' inequality is not violated. At the top, the boundary condition is modified by the hysteresis model adopted for the top trijunction, which is defined by the following generally observed criteria:

1. The contact line does not move, i.e., rc(t) is a constant when the contact angle varies between 4A and OR, where OA and OR are respectively the advancing and receding contact angles (Dussan V. 1979).

2. When the contact angle OA or OR is reached, contact line motion ensues. To simplify the situation, we assume that the contact angle does not vary with contact line motion, i.e.,Oc=OA or kR depending on the direction of motion.

3. It is imposed that an advancing contact line cannot recede without passing through the hysteresis range OA to OR. Thus, if the oscillatory motion of the plate induces an impending recession of an advancing contact line, our model would perforce pin the radius r,(t) at such a location and require 0, to travel from kA to R before allowing recession. Similar conditions hold for the receding angle, OR.

Schematic presentation of the hysteresis model is shown in Figure 3-9. An additional simplification resulting in the model is relevant in the presence of multiple

72

solutions. When the contact line is advancing, for instance, we choose that solution for which r,(t) is closest in the direction of impending motion. So also for a receding contact line. Thus, irrespective of the stability of the equilibrium profiles, proximity of the contact line is made the criterion for choice of a solution from those available. Actually, in the situation where dynamic contact angles exist, energy criteria in terms of the existence of a minimum of the energy curve ( E versus Oc ) is no longer applicable. In particular, the wetting component of energy Ee , namely 7rr c(_YsYsg) cannot be replaced by an expression such as 7rrc',yicosOc, for otherwise ('-y,,yg) would cease to be a material property. Thus, in passing through the successive equilibrium states, the stability of menisci cannot be deduced in terms of an energy minimization principle. From a microscopic viewpoint, if we consider that the hysteresis range 0A to 4R merely represents the apparent contact angle behavior, and that in actuality, the microscopic contact angle is 40 (which is the Young's contact angle, not related in any straightforward way to OA or OR), then we are justified in replacing (,y,--y,,) by -ylcos40. However, on a microscopic scale, the phenomenon of hysteresis is occasioned by the presence of asperities (Oliver et al. 1977, Eick et al. 1975) and surface chemical inhomogeneities (Neumann 1972) or by several other effects which are still a matter of investigation. Thus, the wetted energy expression will need to be modified to -yigcos4OArsj , where Ar., will incorporate the effects of roughness, etc. No tacit assumption regarding E,, may be made either from a microscopic or macroscopic view, and, for this reason, we will refrain from drawing any conclusions regarding stability of the quasi-equilibrium path simulated hereunder.

73

It may be pointed out that the situations presented in this section and the two previous ones represent two limiting types of behavior. In sections 3.5.1 and 3.5.2, the boundary conditions and stability criterion, in terms of the minimization of energy are applicable to a situation where the time scale of motion of the top contact line is greater than that required for the meniscus to achieve equilibrium. Here, on the other hand, it is supposed that the time scale of contact line motion is much shorter than that required for reaching an equilibrium state, the latter being determined only by incorporation of the dynamics of the interface as well as the fluid contained in the meniscus.

Traditionally, the following formula has been employed in models describing growth of a crystal via a pulling process. The variation of radius at the top of the meniscus with pull rate and time is given by (Thomas et al. 1986, Tatarchenko and Brenner 1980),

dr dh
-= (up - S) tan(4- 5) (3.13)
dt d ct

where 0-4, is the deviation in contact angle from the value that leads to fibre growth with non-uniform diameter. We attempt, in what follows, to evaluate the applicability of this expression in the light of hysteresis. In the event of 0 being a constant, it is obvious that the variation of the radius with time follows the same waveform but has an opposite phase with respect to the variation of meniscus height, h,(t). In this respect, the presence of hysteresis introduces a non-linearity in the above expression, in that the value of k in turn depends on r,(t) and its rate of

74

change. However, if a pinning condition arises, i.e., r, = constant, then the above formula indicates that, provided u, ; dh/dt, k = 0., i.e., that the contact angle remains fixed. But this is neither necessary nor true when hysteresis is included, because it can happen that 4. assumes values between OA and R,, while the contact line remains stationary. The results of our calculations, under the restriction of equilibrium, and the conditions 1->3, prescribed by the hysteresis model above, will call to question the implications of Eq. (3.13).

Under the present model, the oscillation of the plate is effected by specifying hc(t)=hc*+Ahsinirt, where in all the results presented, dh,= 0.2h, was given. The hysteresis range 0k to 4R is also specified, with the values corresponding roughly to the material properties of sapphire. We start with an initial condition at the top trijunction which facilitates the calculation of top radius r,(t), when the contact angle is OC=OA+A4A, i.e.,the contact angle is assumed to have been displaced so that the contact line motion may ensue. But for the initial transients, the final periodic behavior of the meniscus was found to be independent of the specified initial state. The subsequent development is followed by plotting the time variation of meniscus profiles, top radius and contact angle. As in the previous sections, Bond numbers of 10' and 2.5x103 are used, and the meniscus is internally pressurized to a given ratio of the atmospheric pressure. Cases simulated below relate to downward pulling of the fibre.

First, a low Bond number case (Bo= 10-) is considered, without hysteresis, the contact angle being fixed at 90*. In Figure 3-10, the sinusoidal motion of the plate

75

given by h,(t), causes a corresponding sinusoidal motion of the top radius, re(t). The r. variation is here out of phase with hc(t) as represented by Eq. (3.13). It is noted that the radius variation in Figure 3-10 for the imposed amplitude Ah, is quite wide, assuming values from 0.62rb to O.83rb. Introduction of a hysteresis range gives rise, in Figure 3-11, to a square waveform for r,(t). For the hysteresis range imposed, OR=850 to OA=950, the pinning condition comes into effect when 0, lies in the hysteresis range. This can be seen from Figs. 3-11 (c) and (d), in which, while 0, varies between 4A and OR, the radius remains a constant. It is not difficult to foresee that when the hysteresis range is further increased, for this amplitude of oscillation, the pinning condition is eventually enforced throughout the cycle of h,(t), while Oc(t) assumes a sinusoidal form within the range OR to 4A In terms of the motion of the contact line, the variation of rc(t) displayed in Figure 3-11 can be interpreted as a stick-slip motion (Jansons 1986).

When we reach a higher Bond number, say 2.5x103 as in Figure 3-12, where the contact angle is fixed at 1000 with no hysteresis and for mean height h, = 0.5, and Ap = 5%. Here the behavior of the meniscus registers a marked change. This is consistent with the difference in meniscus behavior arising due to the different Bond number regimes as seen in the previous sections. It is particularly noticed here that for the same relative variation in height Ah, as in Figure 3-10, in this case the variation in radius rc is less significant. Also the radius variation is in fact in phase with the h, variation, being significantly in opposition to Figure 3-10, and formula, Eq. (3.13). In other words, the radius in this case increases with height instead of

76

decreasing, while the angle remains constant. This is evidently due to the change in meniscus curvature from concave, as in Figure 3-10, to convex. For a higher meniscus, hc* = 1.5, for the same 40 and Ap shown in Figure 3-13, the radius variation is also not very substantial in magnitude and the waveform is not identical to hc(t), indicating the non-linearity that is introduced by 4(t). When hysteresis is present, and for the case of fairly large hysteresis, and tall meniscus, Figure 3-14 shows an interesting situation arising due to the presence of multiple solutions. Note that in the previous cases, only one solution existed in each case, while here 2 or 3 solutions are obtained for certain heights. The choice of solutions is made as explained above. The marked excursions of radius evident in the figure arise when the outer solutions seen in Figure 3-14(a), are allowed to be attained by the imposed conditions 1->3. The hysteresis model in this case causes a noticeable change in the contact angle behavior (Figure 3-14 (d)) in comparison to the static model. There appears to be a Bond number dependence of the phase relationship between dre/dt and dh,/dt, which is not reflected in the relation, Eq. (3.13).

Thus, Eq. (3.13) relating the radius with height of the meniscus, seen in the light of a quasi-equilibrium path, does not appear to exhibit the appropriate features under contact angle hysteresis. In the absence of hysteresis, the Bond number dependence of the phase is not captured. When the Laplace-Young equation yields multiple solutions, the choice of profiles is not straightforward, and the particular selection procedure adopted indicates an abrupt dynamic behavior which may not be justifiably represented by a quasi-equilibrium situation, and the full fluid-dynamic

77

consideration is required to assess the implication of such behavior. We also reiterate that no account has been taken of the stability of the computed equilibrium profiles.

3.6 Conclusions

1. For the boundary conditions adopted here, which pertain some crystal growth processes, namely a fixed base radius, and Neumann boundary condition at the top, solutions to the Laplace-Young equation were obtained using a fourth-order numerical integration scheme. Multiple solutions were obtained for a given contact angle at the top.

2. The energies of the menisci were computed and stability of solutions decided based on an energy minimization concept. While in many cases the energy versus 4c curves did indeed yield a minimum at a specified equilibrium contact angle 00, in other cases local maxima or non-extrema were obtained. Depending on the operating history of the pulling process, it may transpire that an unstable meniscus may form first, as for instance, when the trijunction point of an unstable shape at a new hc/rb is closer to that of the previous equilibrium position.

3. Unless the quantity (h/rb)2Bo is large enough, hydrostatic pressure does not considerably influence meniscus shape, and the effect of pulling direction is therefore not significant. The effect of external pressurization. Ap, also depends on the Bond number. For higher Bond numbers, pressurization can markedly widen the range of

78

available contact angles. Also, for sufficiently tall menisci, pressure (Ap) and gravity (Apgz) interact to cause menisci to assume a variety of curvatures.

4. For modest values of aspect ratio h/rb, there always exist one or more stable meniscus shapes, and Young's contact condition is compatible with the requirement of energy minimization. For higher values of he/rb, although equilibrium solution satisfying Laplace-Young equation do exist, they may not be stable solutions.

5. When the Young's contact condition does not yield a minimum in energy content, then the subsequent non-equilibrium dynamics of the meniscus needs to be simulated. Such treatments have been reported in connection with other situations (Young and Davis 1987, Haley and Miskis 1981) but have not been attempted for the pulling process.

6. Quasi-equilibrium studies of the dynamics of menisci have been made. A simplified hysteresis model was employed to relate the contact angle at the top to the contact line motion. The conditions enjoined by the model enforce pinning of the contact line when the contact angle 0, lies between OR and OA, the advancing and receding contact angles. It is found that the relationship in terms of waveform and phase of r,(t) with respect to hc(t) cannot be captured by the formula usually employed in crystal growth simulations.

79

Flat plate

h,

Meniscus

rb

Pinned at edge

Die

't~c

I

f~ "N 7'/V
/ /

z

f.

Figure 3-1 Schematic of meniscus corresponding to edge-defined fibre growth
process.

r
1

12

80

120 100 80 e 60
40 20

0

(a)

0 20 40 60 80 100 120 140 160 180

(b

10 8 re.rb 6
4

2

V

0.5
0.45
0.4
0.35
0.3
z/rb 0.25
0.2
0.15
0.1
0.05
0

0

I)

20 40 60 80 100

(C)

[K;~

0.5 1 1.5 2

2.5 3 3.5
f/r,

120

4 4.5 5

Characteristics of meniscus formation for Bo = 4x10^3 , AP = 0, hc/rb = 0.5,and upward puiling80
(a) variation of O with angle 4, .
(b) permissible range of r/r, with respect to ,.
(c) profile shapes for different angles O.
(d) profiles corresponding to 4. = 60' .
(e) E profile v.s. 4 , for the specified value of 4 = 60* .
(f) profiles corresponding to 4, = 90* .
(g) E profile v.s. 40 , for the specified value of , = 90*

Figure 3-2

81

2 -1.5 -1 -0.5 0 0.5 1 1.5 2

f/
A0

4.5
4.
3.5
3,
Z-5.
2.
1.5

0.5
0.

0.5
0.45 0.4 0.35
0.3
7/rb 0.25
0.2 0.15
0.1
0.05
0

3
2.8
2-6
2.4
2.2
E 2
1.8 1.6
1.4
1.2

0 20 40 60 80 100 120 140

f)
(i)

(ii)

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2

f'rb

40

(9)

(4H) non-exI2MUM (minimum

0M 40 60 80 100 120 140

Figure 3-2 -- continued

0.5
0.45 0.4 0.35
0.3
0.25
0.2 0.15
0. 1 0.05
0

(d)

(I)

5

E

(e)

(ji) non-exurmum (imminmu

-

82

(a) / (i)

(ii)

10 -8 -6 -4 -2 0 2 4 6 8 10

f/r

410

(b)

)(ii) mamum

2 (C)
1.8 \ (C)
1.6
1.4 (H)
1.2/ /
11 -5 (i)
0.8 0.6
0.4 0.2
0
-5 -4 -3 -2 -1 0 1 2 3 4 5

f/rb

2
1.8 1.6 1.4 1.2 z/rb 1
0.8 0.6
0.4 0.2
0

2.6
2.4 2.2
2

E 1.8 1.6
1.4 1.2

0.8

10 12 14 16 18 20 22 24 26

4d0

10 12 14 16 18 20 22 24 26 28 30

Figure 3-3 Characteristics of meniscus formation for Bo=2.5x106 . AP=5%,
h,/rb =2, and upward pulling.

(a) profiles corresponding to 4, = 16' .
(b) E profile v.s. 6 , for the specified value
(c) profiles corresponding to 4 = 20* .
(d) E profile v.s.4, , for the specified value

of 4o = l6* of o, = 200

5
4.5
4
3.5
E 3
2.5
2
1.5

(d)

(ii) non-CxU-emUm

(i) minimum

83

2
1.8 1.6
1.4 1.2

z~~0.8
0.6
0.4 0.2
0

9.6 9.5
9.4 9.3 E 9.2
9.1
9
8.9
8.8 8.7

40

50 60 70 80 90 100 110 120

(a)

(H)
/(ii) \

-2 -1.5 -1 -0.5 0 0.5 1 1.5
f/rb

10

0i) maxim""' (b)

)minimum

(c)

0)i

-2 -1.5 -1 -0.5 0 0.5 1 1.5 2
f/rb

2
1.8 1.6
1.4 1.2 z/rb 1
0.8 0.6
0.4 0 2

0.015

0.014 0.013 E 0.012

0.011 0.01 0.009

Figure 3-4

Characteristics of meniscus formation for Bo=2.5x103 , ..P=5%, hi/r =2, and upward pulling.

(a) profiles corresponding to $. = 60* (b) E profile v.s.4, , for the specified value of$0 = 60*
(c) profiles corresponding to 0. = 900 .
(d) E profile v.s.40 , for the specified value of 4. = 90*

40 50 60 70 80 90 100 110 120
0 i) minimum

(d)
in mamum

(iii non-ex-mum

2

Full Text

PAGE 1

A STUDY OF FREE AND MOVING BOUNDARY PROBLEMS INVOLVING THIN CRYSTAL GROWTH By Shin-Jye Liang A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 1993

PAGE 2

my parents, my wife, Susan, and my son, Kevin

PAGE 3

ACKNOWLEDGEMENTS I owe an immense debt of gratitude to my parents for the support they have always offered and the love that never wavered. Their belief in me has always been a motivating force behind my achievements. And my wife, Susan, has helped me in more ways than she knows, providing constant support and a few lessons on what is and what should be important in life. I thank my advisor. Professor Wei Shyy, for introducing me to this interesting field of research and providing me with his guidance. His patience and insights have greatly contributed to the quality of this work. Thanks are due to the other members of my supervisory committee. Professors Bruce F. Carroll, Chen-Chi Hsu, Ulrich H. Kurzweg and Craig Saltiel, for their comments and encouragement throughout my time at the University of Florida. Some experimental information used in this study was supplied by Dr. D. G. Backman and Dr. D. Y. Wei; their help is greatly appreciated. I wish also to thank all my colleagues in the CFD group with whom I had helpful discussions, in particular. Dr. Renwei Mei and H. S. Udaykumar for many fruitful suggestions. The present work is partially supported by GE Aircraft Engines, by the National Center for Supercomputer Applications at the University of Illinois, and by the National Science Foundation. iii

PAGE 4

TABLE OF CONTENTS Page ACKNOWLEDGEMENTS iii LIST OF TABLES vi LIST OF FIGURES vii NOMENCLATURE xiii ABSTRACT xvi CHAPTERS 1 INTRODUCTION I 1 . 1 Background 1 1.2 Description of Edge-Defined Film-Fed Growth Technique 2 1.3 Scope of the Present Work 5 1.4 Outline 8 2 INTERACTION OF TIME STEPPING AND CONVECTION SCHEMES FOR UNSTEADY FLOW SIMULATION 13 2.1 Introduction 13 2.2 Model Problem Analysis 16 2.3 Results and Discussions 22 2.3.1 Linear Burgers Equation 23 2.3.2 Nonlinear Burgers Equation 25 2.4 Von Neumann and Spectral Analyses 30 2.5. Summary 33 3 MENISCUS CHARACTERISTICS WITH APPLICATION TO THE EFG PROCESS 53 3.1 Introduction 53 3.2 Configuration, Assumptions and Formulation 56 iv

PAGE 5

3.3 Energy Formulation 60 3.4 Numerics 62 3.5 Results and Discussion 62 3.5.1 Basic Meniscus Behavior with Upward Pulling 63 3.5.2 Effect of Pull Direction on Meniscus Behavior 67 3.5.3 Quasi-equilibrium Motion with Hysteresis 70 3.6 Conclusions 77 4 NUMERICAL TECHNIQUES FOR SOLVING FREE AND MOVING BOUNDARY PROBLEMS WITH APPLICATION TO THE EFG PROCESS 95 4. 1 Introduction 95 4.2 Literature Review 97 4.2.1 Exact Solutions 97 4.2.2 Approximation Solutions 97 4.2.3 Numerical Solutions 98 4.3 Front Tracking and Moving Grid 103 4.3.1 Assessment of the Quasi-stationary Approximation 104 4.3.2 A General Procedure for Interface Tracking 109 4.3.3 Assessment of Time Stepping Schemes for Phase Change Problems Ill 4.4 EFG Process Physics 113 4.5 Mathematical Model 115 4.6 Numerical Simulation of the EFG Process 124 4.7 Conclusions 130 5 DYNAMIC SIMULATION OF THE EFG PROCESS 146 5.1 Introduction 146 5.2 Dynamic Perturbation on the EFG Process with St =1 148 5.2.1 Base Temperature Perturbation 148 5.2.2 Pull Speed Perturbation 151 5.3 Dynamic Perturbation on the EFG Process with St = 0.024 155 5.4 Dynamic Contact Angle 157 5.5 Summary and Concluding Remarks 160 6 CONCLUSIONS AND RECOMMENDATIONS 178 6. 1 Achievements and Findings 178 6.2 Directions for Future Research 181 REFERENCES 183 BIOGRAPHICAL SKETCH ' 193 V

PAGE 6

LIST OF TABLES Tables Page 2-1 Amplification Factors of Backward Euler Scheme for Linear Burgers Equation 35 2-2 Amplification Factors of Crank-Nicolson Scheme for Linear Burgers Equation 36 4-1 Assessment of time stepping scheme for phase change problem. (a) Comparisons of the predicted interface location at r = 0.2 using backward Euler and Crank-Nicolson time stepping schemes for Eqs. (4.3) and (4.4) with 11 and 21 grids. (b) Comparisons of the required CPU time on Micro Vax 3100 using backward Euler and Crank-Nicolson time stepping schemes for Eqs. (4.3) and (4.4) with 11 and 21 grids 132 4-2 Thermophysical properties and physical dimensions of sapphire (AI2O3) used in the present EFG simulations (Case 1: St = 1 and Case 2: St = 0.024) 133 4-3 Definition and magnitude of key dimensionless parameters based on values given in Table 4-1 (Case 1: St = 1, Case 2: St = 0.024) 134 4-4 Two different scaling procedure and resulting nondimensional energy equations and boundary conditions 135 vi

PAGE 7

LIST OF FIGURES Fi gures Page 1-1 Schematic illustration of the EFG sapphire fibre growth equipment used to produce 25 fibre with diameters in the range 0.025 ~ 0. 15 mm. The inset shows a magnified view of the die tip and meniscus region where fibre solidification occurs 11 12 Schematic illustration of the heat transfer mechanisms in the EFG process 12 21 Comparison of solutions of unsteady linear Burgers equation using different convection schemes, backward Euler time stepping, P = 10, C = 1, and 21 nodal points 37 2-2 Comparison of solutions of unsteady linear Burgers equation obtained using the backward Euler and Crank-Nicolson time stepping schemes 38 2-3 Comparison of the backward Euler and Crank-Nicolson time stepping schemes with four convection schemes for unsteady linear Burgers equation using 81 nodal points, and P = 10, C = 1 .... 39 2-4 Effect of time stepping schemes on solution accuracy of the nonlinear Burgers equation with y = 5x10'^ , Ax = 0.02, and At = 0.01 40 2-5 The L-2 error norm with respect to Ax of the nonlinear Burgers equation with p = 5x10"^ and At/ Ax =0.5 41 2-6 Effect of time stepping schemes on solution accuracy of the nonlinear Burgers equation with u = 5x10'^ , Ax = 0.02. and At = 0.01 42 2-7 The L-2 error norm with respect to Ax of the nonlinear Burgers equation with i> = 5x10'^ and At/ Ax =0.5 43 vii

PAGE 8

2-8 Effect of time stepping schemes on solution accuracy of the nonlinear Burgers equation with v = 5x10 ' , Ax = 0.02, and At = 0.05 44 2-9 Performance characteristics of the combined backward Euler and first-order upwind schemes with periodic boundary condition, p = 10* , C = 0.5, and Ax = 0.05 45 2-10 Performance characteristics of the combined backward Euler and second-order central difference schemes with periodic boundary condition, P = 10^ , C = 0.5, and Ax = 0.05 46 2-11 Performance characteristics of the combined backward Euler and second-order upwind schemes with periodic boundary condition, P = 10' , C = 0.5,and Ax = 0.05 47 2-12 Performance characteristics of the combined backward Euler and second-order QUICK schemes with periodic boundary condition, P = 10* , C = 0.5, and Ax = 0.05 48 2-13 Performance characteristics of the combined Crank-Nicolson and first-order upwind schemes with periodic boundary condition, P = 10\ C = 0.5, and Ax = 0.05 49 2-14 Performance characteristics of the combined Crank-Nicolson and second-order central difference schemes with periodic boundary condition, P = 10\ C = 0.5, and Ax = 0.05 50 2-15 Performance characteristics of the combined Crank-Nicolson and second-order upwind schemes with periodic boundary condition, P = 10\ C = 0.5, and Ax = 0.05 51 216 Performance characteristics of the combined Crank-Nicolson and second-order QUICK schemes with periodic boundary condition, P = 10\ C = 0.5, and Ax = 0.05 52 31 Schematic of meniscus corresponding to edge-defined film-fed growth process 79 3-2 Characteristics of meniscus formation for Bo = 4x10 * , AP = 0, hj/rb = 0.5, and upward pulling 80 (a) variation of <^c with angle 0b . (b) permissible range of rjTf, with respect to <^<, . (c) profile shapes for different angles 0b . viii

PAGE 9

(d) profiles corresponding to = 60Â° . (e) E profile v.s.4>^ ,foT the specified value of 0o = 60Â° . (f) profiles corresponding to o = 90Â° . (g) E profile v.s.0,. , for the specified value of 4>Â„ = 90Â° . 3-3 Characteristics of meniscus formation for Bo=2.5xlO'* , AP=5%, hj/rb = 2, and upward pulling 82 (a) profiles corresponding to 0o = 16Â° . (b) E profile v.s.^^ , for the specified value of <^<, = 16Â° . (c) profiles corresponding to (/)<, = 20Â° . (d) E profile v.s.cf)^ , for the specified value of = 20Â° . 3-4 Characteristics of meniscus formation for Bo=2. 5x10"^ , AP=5%, hj/rb = 2, and upward pulling 83 (a) profiles corresponding to o = 60Â° . (b) E profile v.s.<^<. , for the specified value of 0,, = 60Â° . (c) profiles corresponding to 4>Â„ = 90Â° . (d) E profile v.s. 0, , for the specified value of Â„ = 130Â° . 3-6 Effect of the direction of pulling on the meniscus characteristics with Bo=10-^ and A? = 0 85 (a) permissible range of contact angle. (b) permissible range of trijunction location. (c) permissible meniscus shapes at 0,, = 10Â° . (Note that with given Bond number and hJT^, less than 1 , the meniscus shape is dominated by the geometrical constraint from the boundary conditions. Hence the meniscus shapes are essentially the same between two pulling directions. The effect of pulling direction becomes noticeable only when hjTy, becomes larger.) 3-7 Effect of pressurization on the meniscus shapes with Bo =10'^ and A? = 1%. (excess interior pressure) 86 (a) permissible range of contact angle. (b) permissible range of trijunction location. (c) permissible meniscus shapes at 0o = 70Â° . (Note that the hydrostatic pressure effect becomes noticeable only for a higher h,/rb , as evidenced by difference caused by pulling directions.) ix

PAGE 10

3-8 Energy profiles for Bo = 10'^ , AP = 1%, 0o = 70Â° and various hc/rfc. Both upward and downward pulling cases are shown .... 87 3-9 Schematic presentation of hysteresis model 88 3-10 Meniscus behavior at lower Bond number, with no hysteresis. ( Bo = IQ-^ AP = 5%. h/ = 0.5.6. = 6. = 100Â°) 89 (a) Meniscus profiles obtained for the duration of oscillation, (b) Variation of height h^ with time. (c) Computed value of top radius versus time. 3-11 Meniscus behavior at lower Bond number, with small hysteresis range. ( Bo = lO"^ . AP = 5%. h/ = 0.5.6, = 95Â°. 6^ = 85Â°') . . 90 (a) Meniscus profiles obtained for the duration of oscillation. (b) Variation of height h^ with time. (c) Computed value of top radius r^ versus time. (d) Computed variation of contact angle 6c versus time. 3-12 Meniscus behavior at higher Bond number, with no hysteresis. ( Bo = 2.5x10 ' . AP = 5%. h/ = 0.5.6, = 6 p = 100Â°) 91 (a) Meniscus profiles obtained for the duration of oscillation. (b) Variation of height h^ with time. (c) Computed value of top radius r, versus time. 3-13 Meniscus behavior at higher Bond number, with no hysteresis. ( Bo = 2.5x10' . AP = 5%. h-' = \.5.6. = 6u = IQoS 92 (a) Meniscus profiles obtained for the duration of oscillation. (b) Variation of height h, with time. (c) Computed value of top radius r^ versus time. 3-14 Meniscus behavior at higher Bond number, with hysteresis. ( Bo = 2.5x10 ' . AP = 5%. h/ = 1.5.6^ = 110Â°. 6. = 100Â°) ... 93 (a) Meniscus profiles obtained for the duration of oscillation. (b) Variation of height h,. with time. (c) Computed value of top radius r^ versus time. (d) Computed variation of contact angle 6c versus time. 3-15 Meniscus behavior at higher Bond number, with small hysteresis range. ( Bo = 2.5x10 ' . AP = 5%. h,' =0.5.6 . = 90Â°. 6 ^ = 88Â°') 94 (a) Meniscus profiles obtained for the duration of oscillation. (b) Variation of height h, with time. (c) Computed value of top radius r, versus time. (d) Computed variation of contact angle 6c versus time. X

PAGE 11

4-1 Exact solution for melting of a solid confined to a semi-infinite region 137 4-2 Comparison of interface trajectories predicted by the coupled (full equations) and decoupled (quasi-stationary equations) approaches with three values of the Stefan numbers 138 4-3 Comparison of computed and exact solutions for St = 0.1303, St = 1.2216 and St = 2.8576. The figures on the left are the exact and computed temperature fields at different time instants. The figures on the right show the superposed exact and computed interface locations versus time 139 4-4 Schematic configurations and boundary conditions of the present EFG process 140 4-5 Comparison of linear meniscus profile with the LaplaceYoung solution with Rb = 1,H, = 0.188,Rc = 0.763, and various values of Bond number, (a) Ap = 5 % and (b) Ap = 0 % 141 4-6 The grid system and isothermal contours of the initial conditions for both Stefan numbers 142 4-7 The grid system and isothermal contours of the steady-state solution for Case 1: St = 1 143 4-8 The grid system and isothermal contours of the steady-state solution for Case 2: St = 0.024 144 49 Effect of scaling procedure on the convergence characteristics for Case 2: St = 0.024, choice 1: Vn = 0(St) and choice 2: Vn = 0(1) 145 51 Time histories of H,. and subject to a single harmonic perturbation of with (a) Q = 500 At and (b) fi = 20 At (St = 1.0) 163 5-2 Time histories of He and subject to a three-harmonic perturbation of with fi, = 500 At, U2 = 292 At, and ^3 = 108 At (St = 1.0) 164 5-3 The interface shapes at different time instants of (a) a singleharmonic perturbation in ^b with Q, = 500 At and (b) a threeharmonic perturbation ' in ^b with fl, = 500 At, fij = 292 At, and fi3 = 108 At (St = 1.0) 165 xi

PAGE 12

5-4 Characteristics of interface movement with respect to different forcing frequencies of 6^ (a) The relative phase angles between Hj and 6^,, and between and . (b) Percentage variations of H, and R, (St = 1.0) 166 5-5 Time histories and R. subject to a single harmonic perturbation of Up with (a) fi = 500 At and (b) Â« = 20 At (St = 1.0) 167 5-6 Characteristics of interface movement with respect to different forcing frequencies of Up (a) The relative phase angles between and Up, and between H, and R^, (b) Percentage variations of and R, (St = 1.0) 168 5-7 Time histories H, and R^ subject to a three-harmonic perturbation of Up with fi, = 500 At, = 292 At and fij = 108 At (St = l.Q)69 5-8 Time histories of saw-toothed forcing with Q = 20 At and A =0.1, and the corresponding He and R^ (St = 1.0) 170 (a) base temperature perturbation d^ir), (b) pull speed perturbation Up(T). 5-9 Time histories of and R,. subject to a single harmonic perturbation of Up with (a) period =500 At and (b) period =20 At (St = 0.024) 171 5-10 Characteristics of interface movement with respect to different forcing periods of Up (St = 0.024) 172 5-11 A typical Bode diagram depicting sensitivity of the fibre diameter variation with respect to different forcing frequencies of Up (St = 0.024) 173 5-12 Model of dynamic contact condition at the trijunction point (St = 0.024) 174 5-13 Time histories of and R^ subject to a single harmonic perturbation of Up with period = IOOAt (St = 0.024) 175 5-14 Time histories of dH^/dr and dR^/dr subject to a single harmonic perturbation of Up with period = IOOAt (St = 0.024) 176 5-15 Characteristics of interface movement with respect to different forcing periods of Up for different static contact angles (St = 0.024) 177 xii

PAGE 13

NOMENCLATURE Aj amplitude of external forcing of the i* harmonics a capillary length [cm] A-B/A-M Adams-Bashforth/Adams-Moulton predictor-corrector scheme B. E, Backward Euler time stepping scheme Hi Biot number t , Bo Bond number C Courant number C. -N. Crank-Nicolson time stepping scheme Cp specific heat at constant pressure [J/g-K] E total free energy [erg] El gravitational potential energy [erg] E2 surface energy of meniscus [erg] Ej surface energy of the wetted area of the top plate [erg] &2 L-2 error norm EFG process Edge-Defined Film-fed Growth process F 1) nondimensional f(z,t) 2) dimensionless function defining interface x = s(y,t) f meniscus shape distribution as function of r and t [cm] f(x,(^) right hand side of Eq. (2.16) Ge amplification magnitude G, heat flux along the top of the solid phase [W/cm-] g gravitational acceleration [cm/s*] He nondimensional h^ h 1) height of the melt/solid interface h(r,t) [cm] 2) convective heat transfer coefficient [W/cm'-K] he height of the trij unction point [cm] i V-l ij integers denoting grid point in space coordinate Ka wavenumber of the m-th component [cm '] k thermal conductivity [W/cm-K] L nondimensional 1 wavelength of the m-th component [cm] 1 dimensional vertical length of the physical domain [cm] If latent heat of fusion [j/gj Ir reference length scale [cm] Ma Marangoni number N nondimensional unit nbrmal vector xiii

PAGE 14

n unit normal vector along the boundary Nr total number of grid points along the rdirection Nz total number of grid points along the zdirection p 1) cell Peclet number 2) pressure [dyne/cm-] Pe Peclet number Rij Q2> 03. J metrics coefficients defined in Eq. (4-44) R nondimensional radius r radius [cm] Ra Rayleigh number Rad Radiation number radius of die [cm] Te radius of the trij unction point [cm] rf roughness factor (r,z,t) physical coordinates (R,Z,t') nondimensional physical coordinates S nondimensional interface shape S(R,Z,r) = Z H(R,r) = 0 s dimensional function defining interface [cm] St Stefan number T temperature [K] t temporal coordinate [s] T, ambient temperature [K] Tb base temperature [K] Te truncation error for approximating convective term Tj truncation error for approximating diffusive term Tb, melting temperature [K] reference time scale [s] TVD total variation diminishing Up nondimensional pull speed (Up)avg nondimensional average pull speed u convective velocity [cm/s] Up dimensional pull speed [cm/s] Vn normal velocity of interface [cm/s] Wu , Wv nondimensional velocity components in zand rdirection [cm/s] X spatial coordinate [cm] (x,y,t) physical coordinates (X,Y,t*) nondimensional physical coordinates z vertical coordinate [cm] Z nondimensional vertical coordinate Greek a thermal diffusivity [cm^/s] 0 thermal expansion coefficient [1/K] xiv

PAGE 15

r nondimensional thermal diffusivity 7 surface tension . [dyne/cm] A infinite small increment Ap density difference between liquid and gas phases [g/cm^] Ap pressure difference between liquid and gas phases [dyne/cm^] Â€ emissivity e 1) phase angle [deg] 2) nondimensional temperature (^b).vg average base temperature X root of Eq. (4-10) fi dynamic viscosity [g/s-cm] u kinematic viscosity [cm^/s] transformed coordinates p density [g/cm^] a Stefan-Boltzman constant [5.670xlO'^W/cm--K'*] 1) unknown in Burgers equation 2) instantaneous contact angle of the meniscus at the trij unction point [deg] advancing contact angle [deg] receding contact angle [deg] shooting angle at the die tip [deg] contact angle [deg] contact angle at steady-state [deg] app apparent contact angle [deg] Qperiod of external forcing of the i* harmonics [s] superscript * nondimensional quantity subscripts a ambient quantity g gas phase 1 liquid phase m melting point n-1 previous time step n current time step n-f-1 next time step r radius component s solid phase z vertical component XV

PAGE 16

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 A STUDY OF FREE AND MOVING BOUNDARY PROBLEMS INVOLVING THIN CRYSTAL GROWTH by Shin-Jye Liang May 1993 Chairperson: Wei Shyy Major Department: Aerospace Engineering, Mechanics, and Engineering Science Thin sapphire fibers, which find applications in high temperature composites and optical sensors, can be grown from the melt by the Edge-defined Film-fed Gro\\Th (EFG) technique. In this thesis, a thermocapillary model is developed to simulate and illustrate the crystal solidification characteristics. The model, which is axisymmetric and appropriately scaled, solves the energy equation in both the solid and liquid phases separately, and tracks the movement and evolution of their interface explicitly using a newly developed method based on a combined Lagrangian/'Eulerian scheme. The formation of meniscus profiles is examined utilizing the concept of minimization of free energy. Meniscus behavior influenced by Bond number. xvi

PAGE 17

pressurization. aspect ratio, contact angle, and pulling direction has been investigated using both static and quasi-dynamic models for the contact condition. In the course of dynamic simulation, the characteristics of the EFG system subject to base temperature perturbation and pull speed perturbation are studied. To both types of perturbation, the solid/melt interface responds at the forcing frequency. Due to the appearance of multiple time scales in the solidification process, both the model prediction and the experimental observation indicate that as the forcing frequency increases, its effects on crystal growth are reduced. Furthermore, the interface behaviors yielded by the static and dynamic models for the contact conditions at the trij unction point are different, indicating the inadequacy of the conventional approach. The capability developed and insight gained from this work can help improve the EFG process. xvii

PAGE 18

CHAPTER 1 INTRODUCTION 1.1 Background . Free and moving boundaries are encountered in a variety of engineering applications. Analysis of such problems is often difficult because the free and moving boundaries, whose locations are unknown a priori and must be determined as a part of the solution, make the problem nonlinear. The interaction between phase change and other thermofluid transport processes adds additional complexities, making such problems intrinsically difficult to solve. For example, in crystal growth from the melt, the coupled heat transfer and fluid dynamic equations, conjugated with phase change, soUd phase heat conduction and radiative heat exchange with the ambient surfaces, must be solved simultaneously in order to understand the characteristics of the process, the behavior of the growing crystal and the effects of other manufacturing process parameters that determine the quality of the crystal (Tseng et al. 1992, Brown 1988, Viskanta 1988). Mathematically, the free and moving boundary problems are challenging and, except for a few simple cases, have to be solved numerically. Considerable progress has been made in the solution of problems involving free and moving surfaces using numerical methods., A central problem in the numerical modeling 1

PAGE 19

2 of such processes is the treatment of heat and mass flow conditions that arise at the moving solid/liquid interface and the liquid/gas meniscus whose shapes and locations are unknown and vary in time. Many approximation techniques and numerical methods have been developed, as reviewed by Floryan and Rasmussen (1989), Crank (1984), Ozisik (1980), Wilson et al. (1978), Bankoff (1964), Carslaw and Jaeger (1959). The characteristics of the crystal/melt interface is dictated by the Stefan number since it determines the growth rate of the crystal into the melt. In order to control the shape/quality of the crystal, the combined effects of fluid (melt) motion and heat transfer in a variable domain bounded by capillary (melt/gas) and Stefan-type (solid/melt) interfaces need to be studied. Open research problems are of both mathematical (existence, uniqueness of solutions) and practical natures. The main purpose of this study is to numerically analyze free and moving boundary problems involving phase change with special emphasis on the edge-defmed film-fed growth (EFG) process for producing thin sapphire crystals. 1.2 Description of Edge-defined Film-fed Growth Technique Sapphire fibers for use in both optical sensors and structural composites have been produced using the edge-defined film-fed growth (EFG) process (Kalejs et al. 1983a,b, LaBelle 1980, Sachs 1980, Swartz et al. 1975), wherein a crystallographically oriented fibre is grown from a meniscus of molten alumina.

PAGE 20

3 An inert die is used to control the extent of the melt and to help shape the melt and thus the resulting crystal. The EFG process can yield continuous lengths of sapphire fibre with diameters in the range of 0.05 to 0.25 mm. The fibre grower, schematically shown in Figure 1-1 (Backman 1992), contains three main components: the hot zone, the puller system, and the fibre spooling system. The hot zone consists of an inductively heated refractory metal crucible, with the edgedefining dies all enclosed within a water-cooled, environmentally regulated chamber. The melt is supplied to each die tip from a capillary-fed manifold. The puller system consists of a belt puller and a fibre guide. The belt puller utilizes a double belt traction mechanism driven by a precision stepper motor and is monitored with an optical encoder. After the fibre exits the puller, it passes through several pulleys and is wound on spools under regulated tension. Within the hot zone, the fibre growth process begins with seeding at the die tip, which establishes a crystallographic orientation. The puller speed is generally kept constant; the meniscus dimensions and fibre diameter are manipulated through changes in the induction coil power level setting to adjust die tip temperature. In practice, both the die tip temperature and the puller speed may experience either intentional or unintentional variations in time. The unintentional variations can result from the perturbations to the operating conditions caused by fluctuations in the environment, coil power fluctuations, or motor speed irregularities. The intentional variations can be designed to compensate some of the aforementioned perturbations to achieve a controlled growth process and uniform crystal properties and dimensions.

PAGE 21

4 Throughout the development of the EFG process, several problems have been encountered that have potentially severe effects on the quality of the grown crystal. For example, the shape/quality of the thin crystals grown from the melt is strongly influenced by the intricate coupling of heat transfer and melt flow in the EFG system. The heat transport environment must yield a well-controlled temperature field in both the solid and liquid phases so that the crystallization rate, the shape of the solidification interface, and thermoelastic stresses in the crystal can be controlled. A general goal of the design of heat transfer systems for any crystal growth configuration is to establish a nearly constant temperature gradient along the crystal/melt interface and to control the cooling rate of the crystal. Low dislocation and defect densities occur when the temperature gradients in the crystal are low (Kalejs et al. 1983a,b, Sachs 1980 and Swartz et al. 1975). However, because of the small dimensions of the thin crystal sheets, which are of the order of a few mm, the rate of heat losses to the surroundings along the surface of the crystal/melt could be large. This usually results in large temperature gradients near the solidification boundary, which are undesirable. To overcome this problem, the shape-defining dies can be enclosed within a watercooled, environmentally regulated chamber, as illustrated in Figure 1-1. Active control may be desirable to produce crystals of uniform cross section because the shape of the crystal is constrained only by capillary action. A schematic illustration of the heat transfer mechanisms in the EFG process is shown in Figure 1-2. A film of melt forms from the top of a die and is

PAGE 22

5 solidified by convection and radiative heat transfer to the surroundings. It is important to control the growth process according to the heat transfer characteristics in the system to obtain a crystal of definite diameter and desired homogeneity. In the melt, convection driven by mechanisms such as pulling of the seed crystal, buoyancy and surface tension will affect the heat transfer process and hence the solidified crystals. For thin crystal growth from the melt via the EFG technique with a dimension of a few mm, the order of magnitude of convective heat transfer is expected to be relatively small compared with that of the conductive one. Order of magnitude analysis of these heat transfer mechanisms under real operating conditions are evaluated in Chapter 4. Based on this analysis, a mathematical model is developed for simulating the dynamic characteristics of the crystal growth process. Furthermore, the issue of an appropriate scaling procedure, which plays a critical role in determining the computational efficiency, will be addressed. 1.3 Scope of the Present Work The purpose of this study is to mathematically analyze the transport phenomena in the EFG system so as to establish the relationships between the operating conditions and the crystal shape/quality. To develop a complete capability of numerical simulation of free and moving boundary problems involving phase change, several technical issues need to be addressed, including the assessment of time stepping, convection and diffusion schemes, meniscus

PAGE 23

6 behavior and the contact condition at the trij unction location, the solution algorithm for the thermal field in both liquid phase and solid phase, the front tracking technique, and the adaptive gridding to accommodate time-dependent and irregular domain configurations. In order to develop a comprehensive model to conduct dynamic simulations of the EFG process, the following aspects are considered in this thesis. Â• First, the discretization schemes for the temporal and spatial terms and their interaction in the time-dependent thermofluid computations are investigated. In order to clearly understand the issue, both linear and nonlinear Burgers' equations are used as model problems to study the numerical performance of the various schemes. The von Neumann stability and fast Fourier transform techniques are utilized to analyze the stability, accuracy, and dispersive and diffusive characteristics of the combined performance of time stepping and convection schemes. The scope of this aspect of the study is more generic and is intended to provide some general guidelines for helping devise suitable computational strategies for growing crystals of large sizes, where convection effects are more pronounced and may vary spatially. Second, the physics and appropriate treatment of the solid-liquid-gas trijunction point, which is not well understood at the present time, is explored. Consideration of force equilibrium leads to a stress singularity at the contact line. The modeling of static contact lines, where a contact angle is specified at the trijunction point, is justified only in the liquid-gas-solid system under the

PAGE 24

7 assumption that the solid wall is smooth and homogeneous. The extension of this concept to dynamic contact lines is desirable, to relate the instantaneous contact angle to the velocity of contact line motion. The attempts to determine this function from physical experiments have not been very successful and this remains an open area of research. Relatively detailed discussion will be given later, for both the static and dynamic contact conditions. A simplified dynamic model, taking into account the effect of the direction and speed of contact point movement, will also be proposed and incorporated into the present model. Based on both the static and dynamic contact conditions, meniscus characteristics and solidification behavior in the EFG process will be examined and compared. The third aspect of this study is to simultaneously solve the equations of heat transfer in both solid and liquid phases, liquid/solid interface movement, and the dynamics of meniscus and contact lines using a curvilinear coordinate system. Body-fitted coordinates are used for mapping the irregular shape of the timewisevarying solid/liquid and liquid/gas interfaces. A new Lagrangian approach is developed for tracking the phase front, and a solution procedure is used for solving simultaneously the thermal fields in both liquid and solid phases. This yields a complete EFG numerical simulation. Finally, thermal sensitivity and stability of the EFG for sapphire fibre production is numerically studied. The dynamic behavior of the EFG system with respect to an external disturbance, e.g., varying the base temperature or pull speed, is investigated. The proposed dynamic model is incorporated to study the

PAGE 25

8 dynamic EFG characteristics. Based on this model, the variation of radius of the grown crystal is compared with the static wetting condition, which is conventionally employed in numerical simulations of crystal growth (Crowley 1983, Ettouney and Brown 1983). The relationship between response and forcing frequency, in terms of magnitude and phase shift, is identified. 1.4 Outline The present study therefore can be divided into several parts: discretization schemes for temporal and spatial terms and their interaction, meniscus and contact point behaviors, front tracking and solution techniques for phase change problems, and dynamic simulation of the EFG process. In Chapter 2, the time stepping schemes, the convection schemes, and their interaction in unsteady thermofluid computations are presented. Both linear and nonlinear Burgers' equations are used as model problems to study the stability, accuracy, and dispersive and diffusive characteristics of each scheme. With the aid of a stability analysis, the overall performance of the scheme is discussed. The shape of the melt/gas meniscus that connects the growing crystal with the die, which is governed by the LaplaceYoung equation, is critical in determining the properties of the solidified crystal. In Chapter 3, the meniscus characteristics with application to the EFG process are studied. Both static and dynamic contact line motions are investigated. Quasi-equilibrium dynamics of the meniscus is simulated using a simplified hysteresis model for the contact angle at

PAGE 26

9 the trij unction point. The implication of the variation of the meniscus shape on the crystal growth characteristics is explored in detail. In Chapter 4, a general treatment of the interface movement is designed, which is not limited to a single-valued or isothermal surface employed by other investigators (Lacroix 1990, Thompson and Szekely 1989, Ettouney and Brown 1983, Crowley 1983, Sparrow et al. 1977). The development of the numerical technique is reported in Shyyet al. (1992b), which demonstrates that the method is robust and can handle both isothermal and nonisothermal phase change problems. The technique is combined with a moving gridding procedure, determining the thermal field in both the solid and liquid phases and the evolution of both the solid/melt interface and the melt/gas meniscus to predict the solidification characteristics of sapphire fiber growth via the EFG technique for different Stefan numbers. Two different scaling procedures for nondimensionalizing the governing equations are examined, and the computational characteristics and relative efficiency of these two procedures are compared. In Chapter 5, the sensitivity and dynamic characteristics of the EFG system subject to perturbations of the operating environment are explored. The dynamic behavior of the EFG system with respect to the external disturbance, i.e., variations in either base temperature at the die tip or pull speed, is investigated. A dynamic model is also incorporated into the system and compared with the conventional static model for the contact condition at the trijunction point. In Chapter 6, a summary and some suggestions for further work are given.

PAGE 27

10 Overall, the contributions of this thesis are to develop and utilize advanced and original numerical techniques for simulating the EFG process, to examine and assess the contact conditions relevant to the present physical problem, and to construct a predictive tool incorporating both the physical models and the numerical techniques. The developed model system is envisaged as a computational tool, to make detailed investigations of the dynamic characteristics of the EFG processes.

PAGE 28

11 4J U -S j= Â£ <^ ^ c = . o ^ Â£ .2 E 5 w 00 CO x: sj o Â£.Â•1 a. > Â— pre CO 3 :S S I 2 ^ ^ ^ Â« S ^ .1 = :2 s .Si c O I o cr Â• c CO u O CO Oi)

PAGE 29

12 pulling speed convection & radiation Solid (diffusion) interface (latent heat)^ Liquid (diffusion & convection) trij unction point meniscus (Young-Laplace Eq) Figure 1-2 Schematic illustration of heat transfer mechanisms in EFG process.

PAGE 30

CHAPTER 2 INTERACTION OF TIME STEPPING AND CONVECTION SCHEMES FOR UNSTEADY FLOW SIMULATION 2.1 Introduction The numerical solutions of heat transfer, fluid flow, and other related transport problems can be obtained by solving the laws governing these processes in discrete forms. By using appropriate numerical techniques to discretize the spatial derivatives, the original continuous problems are reduced to a system of ordinary differential equations in time. A time stepping scheme can then be used to advance the solution in time from given initial and boundary conditions. Categorically, two kinds of time stepping methods can be defined, namely, explicit and implicit. The explicit schemes provide a noniterative "marching" procedure for obtaining the solution of each nodal point at present time in terms of the known preceding and boundary values. Explicit schemes limit the time step size due to the stability as well as the accuracy requirements. Implicit procedures, on the other hand, generally involve the simultaneous solution of all unknowns in terms of both the known preceding and the unknown present values at other nodes. The implicit methods are usually more stable, and their restriction in time step size is mainly due to the accuracy consideration. 13

PAGE 31

14 Besides the time stepping methods, the convection schemes play an important role in both steady and unsteady flow computations. Discretization of the convection terms is one of the principal difficulties in computing the complex fluid flows. The first-order upwind scheme used to be a popular choice due to its superior numerical stability as well as the wiggles-free characteristics. However, the resulting accuracy yielded by this scheme is often found to be inadequate in view of its first-order truncation error and strong numerical, as opposed to physical, diffusion effects (Shyy 1985,Roache 1972). Many alternative convection schemes have been studied, as reviewed in Hirsch (1990), Fletcher (1988), Correa and Shyy (1987), and Warming et al. (1973), and the findings for the same schemes are not necessarily consistent among different studies; they are obviously problem as well as implementation dependent. Furthermore, it is known that the performance of a given convection scheme can be quite different between the steady and unsteady flow problems (Shyy 1984), hence, it is important to address the issues related to the interactions between the time stepping and the convection schemes for unsteady flow calculations. In this chapter, assessment of the relative merits of the various possible choices for solving the convection-dominated transport problems is made through a series of well defined cases. The viewpoint taken here is that a numerical scheme preferably is of high order of accuracy in truncation error not only so that it can perform well for problems of smooth solution profiles, but, equally importantly, so that no over/undershoots or spurious oscillations that violate the conservation laws appear, in order to

PAGE 32

15 be able to honor physical realizabiiity. In the latter regard, the concept of total variation diminishing (TVD) has been found very helpful for constructing an appropriate convection scheme (Hirsch 1990). It is known that, with the use of the forward Euler scheme, the only convection scheme that can guarantee the TVD property across a sharp gradient or a discontinuity is first-order (Roe 1986). Here, through pragmatic testing and basic analysis, the intent is to gain insights of the performance of the schemes considered without insisting on the satisfaction of TVD. The reasons are twofold. No known technique is available to enforce TVD for a uniformly high order time stepping scheme; furthermore, for problems with substantial source terms, the concept of TVD does not apply (Hirsch 1990). Three time stepping schemes, the first-order backward Euler, the second-order Crank-Nicolson. and the third-order Adams-Bashforth/Adams-Moulton (A-B/A-M) predictor-corrector schemes, are considered. In conjunction, four convection schemes, first-order upwind, second-order central differencing, second-order upwind, and QUICK (Leonard 1979), are examined using both the linear and nonlinear Burgers equations. The goal is to understand the dispersive and dissipative characteristics exhibited by each scheme and to offer guidance for utilizing these features to help improve the overall numerical accuracy. The Burgers equation, of the linear as well as nonlinear forms, is used as the test problem. Furthermore, the von Neumann stability analysis and the FFT spectral analysis will also be used to help understand the characteristics of each formulation.

PAGE 33

16 The results to be presented in the following are a part of the attempt of developing the necessary models for predicting the crystal growth processes with different physical scales and operating conditions. As will become clear later, the convective effect is not dominant in the present EFG processes; however, the work conducted here will become essential for any system of substantial convective effects. Assessment of the time stepping schemes for phase change problem, which is critical for the problem under study, is presented in section 4-3. 2.2 Model Problem Analysis The model problem considered is the unsteady one-dimensional Burgers equation ^^^{u4>)=Â—{i>Â—4>), 00, v=comtant >Q (2.1) at dx dx dx where if u is prescribed constant, the equation is linear, while if u is taken to be a function of (f), then it is nonlinear. Both cases are adopted here since they are relevant to practical applications under different circumstances. The boundary and initial conditions are given in the following : (I) Linear Problem ( u = 1 ) (i) boundary conditions 0(0,0 = 1 (l,t) = 0

PAGE 34

17 (ii) initial condition (x,0) = (II) Nonlinear Problem ( u = <^ ) (i) boundary condition 1 , for 0ix,l) = l+(-)^exp(Â£!) where to = exp(l/8v). The corresponding analytic solution (Lohar and Jain 1981, Benton and Platzman 1972) is X Hx,t) = i (2 IMÂ— )'exp(Â£l) In what follows, the second-order central differencing scheme is always adopted to approximate the diffusion term, i.e..

PAGE 35

(2.3) where T,, , of the order of 0[(Ax)^], is the truncation error. As to the convection term, the following schemes are considered and tested: (a) first-order upwind d(u4>) dx (u)^ (M<^),., Ax Ax + r , foru>0 + r . foru <0 where TÂ„ of order of 0[(Ax)], is the truncation error, (b) second-order central differencing d(u({>) I ^ iu\., ^ dx 2Ax where T, is 0[(Ax)^] (c) second-order upwind (2.4) (2.5) d(u) _. dx -L[ 3(M0), -4(M0),., +(M<^)...3 ] + r , foru>0 ZAX (^.O) 1 2Aa: where T, is 0[(Ax)^]

PAGE 36

19 (d) QUICK diu4>) I ^. dx J_[ -{uci> ) . +7(u<^. ) . -3(M<;.) . -3(M<^), +7 ,/or u < 0 8Ax where is also 0[(Ax)^]. When the first-order backward Euler method for the time derivative is used, the discretized form of Eq. (2.1) can be expressed as follows: (a) first-order upwind where P = u Ax/p = cell Peclet number, C = u At/ Ax = Courant number and the superscript "o" represents the known value. (b) second-order central differencing (c) second-order upwind

PAGE 37

20 (d) QUICK (2.11) The Crank-Nicolson scheme is a second-order accurate method and usually is described as unconditionally stable (Fletcher 1988, Anderson et al. 1984). However, as demonstrated later, this does not mean that the solutions are free of spurious oscillations. If the convection and diffusion terms are replaced by the appropriate discretization schemes at both the current and previous time levels, the corresponding discretized equations result: (a) first-order upwind (2.12) (b) second-order central differencing (Ul) 0,.. -(2-2Z) 0,. .(1-^) 0,.^, = (2.13)

PAGE 38

(c) second-order upwind 21 (2.14) (d) QUICK (2.15) The third-order predictor-corrector scheme, based on the Adams-Bashforth scheme as the predictor step and the Adams-Moulton scheme as the corrector step, is explicit; it involves successive computations at each node but does not need to solve the whole domain simultaneously. If Eq. (2.1) is rewritten in the following manner, at ox dx dx (2.16) then the present predictor-corrector algorithm results : Adams-Bashforth method : (2.17)

PAGE 39

Adams-Moulton method 22 ^c'l) = + ^( 5/c*i)+ 8/w ) (2.18) where (n+1)* designates the value obtained from the predictor step. It is noted that the appropriate convection and diffusion schemes can be used to represent the function f(x,<^). The present A-B/A-M scheme needs a starting procedure since it requires information from the previous two nodes. Because the solutions close to the left boundary are smooth, it is found that no discernible difference exists between using the Crank-Nicolson and the backward Euler schemes as the starting procedure. 2.3 Results and Discussions The viewpoint taken here is that a numerical scheme preferably is of high order of accuracy in truncation error not only so that it can perform well for problems of smooth solution profiles, but equally importantly, so that no over/undershoots or spurious oscillations that violate the conservation laws appear, in order to be able to honor physical realizability. In the latter regard, the concept of TVD has been found very helpful for constructing, as well as assuming, the convection scheme Hirsch (1990). It is known that, with the use of the forward Euler scheme, the only convection scheme that can guarantee the TVD property is first order. Here, through pragmatic testing and basic analysis, the intent is to gain insights of the performance of the schemes considered without insisting on the satisfaction of TVD. The reasons are twofold. No known technique is available to enforce TVD for a

PAGE 40

23 high order time stepping scheme; furthermore, for problems with substantial source terms, the TVD scheme does not work well anyway (Hirsch 1990). The results of the linear and nonlinear Burgers equations are presented in the following sections. 2.3.1 Linear Burgers Equation The linear Burgers equation is chosen because it is relevant to many heat transfer applications. It is known that at steady state, both first-order and secondorder upwind solutions are wiggles-free while neither second-order central differencing nor QUICK is wiggles-free (Shyy 1985). Their relative merits in the context of the unsteady problems are assessed here. The time stepping schemes employed here are forward and backward Euler schemes. The A-B/A-M scheme will be included later in the nonlinear case. For the linear problem, with C = 1 and P = 10, either 21 or 81 nodal points are utilized in the computations. The results with 21 nodes are first presented and compared in Figures 2-1 and 2-2 for the two time stepping schemes at t = 0.3,0.8, and 2. Using the backward Euler scheme, as shown in Figure 2-1, both the firstorder and second-order upwind solutions are wiggles-free, but the first-order upwind is a little more dissipative. For the steady state solutions of Eq. (2.1), it is well known that the central differencing scheme exhibits oscillations if the cell Peclet number, P, exceeds 2 (Fletcher 1988 and Shyy 1985). Figure 2-1 shows that this scheme also produces nonphysical oscillations behind the sharp front for the unsteady

PAGE 41

PAGE 42

25 scheme, Crank-Nicolson, being less dissipative but yielding no wiggles, performs better than backward Euler. Similarly, the performance of a second-order time stepping scheme such as Crank-Nicolson is less satisfactory when combined with other second-order convection schemes due to the fact that both time stepping and convection schemes are dispersive. As illustrated in Figures 2-2 (b) -(d), nonphysical solutions are observed for all three second-order convection schemes. For both the central differencing and QUICK schemes, the oscillations yielded by Crank-Nicolson are substantially larger than those by backward Euler. For the second-order upwind scheme, while the Crank-Nicolson scheme yields improved solutions with sharper gradients and no oscillations at t = 0.8, it produces a noticeable nonphysical undershoot at t = 0.3. The results presented in Figures 2-1 and 2-2 are generic and not affected by the variations in the number of grid points. In order to demonstrate this fact, Figure 2-3 compares the performances between the Crank-Nicolson scheme and the backward Euler scheme, along with the four different convection approximation schemes, where the values of P and C remain the same as before (10 and 1, respectively), but the number of nodal points is increased from 21 to 81. All the major characteristics depicted in Figure 2-3 can be observed in Figure 2-2. 2.3.2 Nonlinear Burners Equation For the nonlinear Burgers equation, two different values of viscosity, v = 5x10-' and 5xl0\ are used to explore the effect of solution gradients on the

PAGE 43

26 performance of each scheme. Three time stepping schemes, including the first-order backward Euler scheme, the second-order Crank-Nicolson scheme, and the thirdorder Adams-Bashforth/Adams-Moulton predictor-corrector scheme, are studied along with the first-order and second-order upwind schemes. Solutions obtained at two time instants, i.e.,t = 2 and 3.5, are presented. The case of v = 5x10"^ is considered first. Figure 2-4 compares the solutions yielded by the three time stepping schemes in conjunction with (a) the first-order upwind convection scheme and (b) the second-order upwind convection scheme. All the computations shown in Figure 2-4 are based on 51 nodes (Ax = 0.02), and At = 0.01, resulting in P = 4u and C = 0.5u. Figure 2-4 shows that with first-order upwind, the numerical solutions yielded by backward Euler and Crank-Nicolson are fairly comparable while noticeable improvement can be obtained by using the thirdorder A-B/A-M scheme. Figure 2-4 indicates that a first-order scheme in space and a third-order scheme in time can actually form a good combination. Regarding the use of the second-order upwind convection scheme, shown in Figure 2-4 (b), all three time stepping schemes perform well, producing solutions in good agreement with the analytical one. In order to measure more quantitatively the performance of each individual scheme under study, Figure 2-5 depicts an L-2 error norm with respect to Ax, where the Courant number is held fixed for all cases, i.e.. At and Ax vary in exact proportion. The L-2 error norm, CjCt), is defined as follows,

PAGE 44

27 where
PAGE 45

28 demonstrate that the combined performance of the discrete operators both in time and in space may not be the same as that suggested by a local truncation error analysis of the individual operator. The case of = 5x10"^ represents a problem of relatively smooth solution profiles. A more stringent case of v = 5x10'^ is considered next. Figure 2-6 compares the numerical solutions corresponding to the different combinations of the convection and time stepping schemes. Comparing Figures 2-4 and 2-6, it is clear that the degrees of satisfaction yielded by these schemes between the two values of u are very different. With a reduced u, the solution gradients stiffen, and the numerical wiggles are more prone to appear. Figure 2-6 (a) shows that, even with the use of the first-order upwind convection scheme, which is viewed as usually too dissipative, nonphysical wiggles develop when the A-B/A-M time stepping scheme is employed. With either backward Euler or Crank-Nicolson, on the other hand, the solutions are of smooth, diffusive characteristics that exhibit no overshoots. The difficulties of resolving the sharp gradients become more pronounced for the secondorder upwind scheme, shown in Figure 2-6 (b), in which the numerical wiggles appear in all solutions. The L-2 error norms of the solutions of u = 5x10'* are summarized in Figure 2-7 where the Courant number again remains unchanged as Ax is varied. Unlike those curves shown in Figure 2-5 for v = 5x10'^, the results in Figure 2-7 for p = 5x10"' are of inconsistent appearances, depicting no recognizable correlations with Ax among different discretization schemes or between different time instants. In

PAGE 46

29 particular, the solutions associated with the A-B/A-M scheme are most erratic; they can either be less accurate than the solutions using other time stepping schemes or produce more error as Ax is reduced. The results presented so far clearly suggest that the three time stepping schemes exhibit different dispersive characteristics; as the leading local truncation error of a time stepping scheme becomes higher-order in At, the resulting solutions tend to be more dispersive, in form of either long wave overshoot or short wave wiggles. The above discussions are based on a single set of the Courant and cell Peclet numbers. Since each of the time stepping and convection schemes possesses its own diffusive and dispersive characteristics, it appears useful to adjust the values of P and C so that a time stepping scheme of particular dispersive and diffusive characteristics can be utilized in conjunction with a convection scheme of complementary characteristics to improve the overall solution accuracy. That such a possibility does exist is illustrated in Figure 2-8, where both 5xl0^and Ax, 0.02, remain the same as those used in Figure 2-6, while At is increased from 0.01 to 0.05, resulting in a C = 2.5u. By increasing the size of At, the backward Euler scheme exhibits stronger diffusion characteristics while the other two higher-order time stepping schemes exhibit stronger dispersive characteristics. Hence, as demonstrated in Figure 2-8 the first-order upwind convection scheme, being generally more diffusive, can be used along with the Crank-Nicolson scheme, a more dispersive scheme. The net outcome is that they can produce satisfactory solutions exhibiting no spurious wiggles with a larger time step, yielding saving in the computing cost. On the other hand, the

PAGE 47

30 second-order upwind convection scheme, being less diffusive, can apparently benefit from the backward Euler time stepping scheme, which is more diffusive. These interesting features suggest that the different characteristics in convection and time stepping schemes can indeed be selectively combined to form an overall algorithm capable of suppressing excessive numerical dispersion and/or diffusion with reasonable At. 2.4 Von Neumann and Spectral Analyses In order to gain some insight of the dissipative and dispersive characteristics of a given numerical scheme, the von Neumann stability analysis (Warming and Hyett 1974, Morton 1971) and the discrete fast Fourier transform (FFT) spectral analysis are utilized. The von Neumann stability analysis is applicable to linear, constant coefficient problems with periodic boundary conditions. The elementary solution for linear Burgers equation is = g (2.20) where i = V-1, kÂ„ = 2ir/LÂ„ is the wave number of the m-th component, and LÂ„ is the wavelength of the m-th component. Based on the elementary solution, the amplification factor of the exact solution can be defined

PAGE 48

31 G = ^(^^^) = e'^^'e-'^"" (2-21) 0(0 where | Gj | = exp(-C/3Vp) is the amplification magnitude, 0 = jSC is the phase angle, jS = kÂ„ Ax, and 27rAx//3 is the wavelength. It is clear that depends on the cell Peclet number, Courant number, and the wave number. The amplification factors yielded by the von Neumann analysis for the backward Euler and Crank-Nicolson time stepping schemes along with various convection schemes for linear Burgers equation are given in Tables 1 and 2. respectively. In order to illustrate the characteristics of each scheme more clearly, a series of results are computed for P = 10^, C = 0.5, Ax = 0.05, and a periodic boundary condition. The initial condition is defined as _/l, /or 0.05<^ <0.2 or 0.55<;c <0.7 '^^'''^^ = {o, themlse Each of Figures 2-9 to 2-16 includes four parts : (i) the comparison of exact and numerical solutions at t = 2At, (ii) the spectral analysis of the numerical solutions yielded by FFT, and the comparison of (iii) the amplification factors as well as (iv) the phase angles according to the von Neumann analysis. The performance of the backward Euler time stepping scheme along with different convection schemes are presented in Figures 2-9 to 2-12; the corresponding results of the Crank-Nicolson time stepping scheme are presented in Figures 2-13 to 2-16.

PAGE 49

32 For the linear problem subject to a periodic boundary condition, the ratios of magnitudes at each wave number between two consecutive time steps shown in the spectral plot (iii), are identical to the amplification magnitude shown in (iii). It is noted that with a high Peclet number, P = 10\ and a periodic boundary condition, the exact solution essentially maintains its initial profile just like a pure wave equation. Consequently, the spectral distributions of the exact solution remains unchanged with time,and the amplification magnitude shown in (iii) are equal to 1 at all wave numbers. Several observations can be made by inspecting Figures 2-9 to 2-16. First, Crank-Nicolson does not yield less dissipation than the backward Euler scheme for all wave numbers. In fact, except for the case of the central difference convection scheme, part (iii) of Figures 2-9 to 2-16 demonstrate that for all three upwind type of convection schemes, Crank-Nicolson is less dissipative than backward Euler only in the long wave length regime; with shorter wave, it is actually more dissipative than backward Euler. Secondly, the amplification magnitudes of the second-order central difference scheme are very close to the exact values; when combined with Crank-Nicolson, it is indistinguishable from the exact profile, i.e. there is no numerical damping created by the scheme, which is good. However, with the central difference scheme, the short waves travel with much smaller phase speeds than the exact solutions. It is these phase speed errors in the short wavelength that cause the central difference scheme to perform unsatisfactorily. This observation explains why the central difference scheme usually yields solutions with highly noticeable 2Ax oscillations.

PAGE 50

33 The other interesting aspect revealed by the present analysis is that the second-order upwind scheme is actually more dissipative than the first-order upwind scheme at short wavelengths as indicated by part (iii) of Figures 2-9 and 2-11. With both Crank-Nicolson and backward Euler, the second-order upwind scheme is less dissipative than the first-order upwind scheme only in the long wave regime. This feature may be useful if a solution profile contains sharp gradients because the phase speed errors associated with short waves are less pronounced. Furthermore, by comparing Figures 2-11 and 2-15, it becomes clear that it is the phase errors in the long wavelength that are mostly responsible for the inaccuracies created by the combined Crank-Nicolson/second-order upwind schemes. This can explain why its solutions exhibit a single long wave overshoot instead of 2Ax short wave oscillations. For the first-order upwind scheme, it appears that in the long wave regime, CrankNicolson is not only less dissipative but also more accurate in phase speed. Accordingly, as already concluded previously, backward Euler/ second-order upwind or Crank-Nicolson/first-order upwind form good combinations. The performance of QUICK lies between the central difference and upwind schemes. 2.5 Summary Both linear and nonlinear Burgers equations have been used as the model problems to investigate the interaction of the time stepping and convection schemes for unsteady flow simulation. The von Neumann stability analysis and the FFT spectral analysis are also utilized to aid our understanding of the performance

PAGE 51

34 characteristics of these schemes. Interesting features are revealed by these analysis. For example, Crank-Nicolson is less dissipative than backward Euler only in the long wave regimes; it is actually more dissipative as the wavelengths are shortened toward 2Ax. Similarly, with the same time stepping scheme, the second-order upwind scheme is less dissipative than the first-order upwind scheme mainly for the longer waves. Consequently, for the second-order upwind scheme, backward Euler yields more accurate phase speed than Crank-Nicolson; while for the first-order upwind scheme, it is Crank-Nicolson that is more satisfactory in phase speed. Regarding the second-order central difference scheme, because it contains very little numerical damping for all wavelengths, the large phase speed errors, especially at short wave lengths, cause its solutions to exhibit pronounced 2Ax short wave oscillations. In contrast, for second-order upwind, because of its higher dissipation rate at short wavelengths, it is the phase errors of longer waves that cause the noticeable inaccuracies; consequently, the solution error tends to appear as a single long wave overshoot. The performance of QUICK essentially lies between the central differencing and upwind scheme. Overall, among all the schemes tested, either a combination of first-order upwind for convection and Crank-Nicolson for time, or a combination of secondorder upwind for convection and backward Euler for time performs better. A consistent second-order accuracy for both time and space discretizations does not produce good results for the time dependent flows with large gradients. The results suggest that a selective combination of the convection and time stepping schemes of appropriate dispersive and diffusive characteristics can attain good performance.

PAGE 52

35 a. E to CO a. a, + A. w I a,|(j .S CO 0. 8 CM + + I a,|0 CO. .S CO a. cs. a,i 00 Q. U E Â— Â«

PAGE 53

36 c _o u E CO a. en. u ST + I .5 V) a. a, + + tN + a. + CM CO. a. If I I CN ca. .5 CO a. I ca. a,Tu + I ca. .S .S 09 a,|cs I ta. + caL CM a,|rj + aTfu CM I a. + .s a, CN I ca. ^ I C5 a,|cN + CN ca. a.|co + a-Tu CN I a, cn |oo + CN ca. .S CO ca. CN .s' a,|co u aTTcN + ea. CN 8 a,|oo a,|0 CN + a, cn I oo + CN I II E u E u. U U f a u k. (/] o CO C Â£ S Q. 3 U u. C T3 C u (J u II r Â« o .5 M D. 3 U E Â— Â«

PAGE 54

37' 5 s II d in o S Â° " T Â« O 01 o <-> S Â•u ^ = i3) -o tn a CM d o ^ Â« II u 5 c o CO 3 Wl 00 d (A O u j= o II c ISO on C '5 w f8 u Q. c u E > "(3 c Â•o o U CO no CM a s on

PAGE 55

38 Â— Backward Euler Crank-Nicolson V E u c 'o c u a Urn c u u u Â•13 13 c o u u u E u u 0. 3 13 1Â« a II ( r> / 0 / 1 T3 00 d u E u u y 5 u c Q. 3 1) "O o I Â•o c o u u CO 00 e Â•8 O 4J Â§ J Â— o 5* 00 Â« c en & ^ Q. C ca *" u = 1 >Â» o o o
PAGE 56

39

PAGE 58

41

PAGE 60

43

PAGE 61

44 4) U > C o u 5 a. c o u u u E u > c o o a. o o O O ^ II L5 < u CO "O c o d II X <3 c o 13 w-i OA ^ Â•|:Â£ u o 5 uu aa oo t u w 3 00

PAGE 62

45

PAGE 63

46

PAGE 64

47 d '^ o d o a IP

PAGE 67

50

PAGE 69

52
PAGE 75

58 increases with Bo, since a higher Bond number implies a taller meniscus. When Ap?iO, the effect of Bond number will influence menisci through the relative values of hydrostatic pressure, Apgz, and externally imposed pressure difference, Ap, as will be shown from our calculations. The Young's contact condition at the top of the meniscus is, where 7,g, 7,, and 7,5 are respectively the surface tensions of the solid-gas, solid-liquid and liquid-gas interface. is the static contact angle measured with respect to the liquid-solid interface. It may be noted that this relation is strictly applicable to a perfectly smooth surface without energy heterogeneity. It merely defines the static contact condition at the three-phase line for an ideal surface. However, if the solid surface is taken to consist of 'patches' of ideal smooth surfaces, we may impose local Young's contact angle at the contact line. This situation has been shown previously to offer one plausible model for hysteresis (Huh and Mason 1977, Eick et al. 1975, Neumann and Good 1972). It is also evident that such heterogeneities or those associated with crystalline orientations or defects, may cause the meniscus to deviate from axisymmetry. Contortions of the contact line are developed in order to ensure energy minimization, by locally enforcing on the liquid-gas surface the Young's angle at the contact line. This in turn distorts the meniscus surface which, for equilibrium, has to conform locally to the Young's condition. In our work, we avoid the complications associated with these aspects by assuming an ideal surface at the top. It is a

PAGE 76

59 conventional practice to account for the contact angle hysteresis by adopting a relation such as Wenzel's, namely where rf is a roughness factor and (i>,Â„ is the apparent contact angle (Oliver et al. 1977). So as not to dilute our focus on the fundamental meniscus behavior under equilibrium conditions, this practice is not adopted. It is possible, however, by scanning through the available range of 4>o for solutions to the LaplaceYoung equation, to incorporate the effect of Eq. (3.6) by resorting to the appropriate value of 4>,^. Hence, a dynamic process with varying (^,pp can be treated simply as a juxtaposition of equilibrium solutions with varying 0<,The LaplaceYoung equation being an ordinary differential equation of second order, two boundary conditions are necessary. For the EFG process, the radius is fixed at the base. The liquid-vapor line is assumed to be pinned at the edge of the die as shown (Surek et al. 1980, Oliver et al. 1977), such pinning having been observed in practice. In the present simplified framework, i.e., without explicitly accounting for heat flow and solidification, we choose to picture the meniscus as attached to a flat plate. Thus, at the top, for a given meniscus height, h, , the contact angle is specified and the radius of the trij unction point can be computed. This is in cognizance of the fact that the Young's contact angle is a material property, independent of the geometry of the meniscus. Such a boundary treatment also offers a simple way to simulate fibre diameter variations. However, the fact that a Neumann boundary condition is

PAGE 77

60 specified at the top results in the existence of multiple solutions. In fact, it is known that a boundary value problem generally does not guarantee a mathematically unique solution. Hence, guidance will be needed to help select the solution that is physically stable. The implications of mathematically non-unique solutions on the realizability of the meniscus shape will be discussed later. 3.3 Energy Formulation The thermodynamics of the meniscus may be studied by defining the energy change in forming the meniscus. The total energy E of this system is given by the sum of the gravitational potential energy and the free energy of the surface. To nondimensionalize energy, Â£Â• = Â— ^ (3.7) ^pga* Total energy resulting from the formation of the meniscus, E, is equal to the potential energy (from the effective head) E, plus surface energy of meniscus E, plus surface energy of the wetted area of the top plate E,. The first component of the energy, that due to gravity and imposed pressure is, Â£i = fitf^ (Apgz+A/?) ek (3-8) The second component, corresponding to that expended in forming the free surface is,

PAGE 78

61 (3.9) The third component, the energy of wetting of top surface is, ^3 = (YfaY,^) T^r] (3.10) where r, is the radial location of the trijunction at the top of the meniscus. We follow the usual practice of replacing the difference, (7,, -7,5) in the third component using the Young's contact condition (Pitts 1976, and Huh and Scriven 1960). This procedure facilitates evaluation of this component, since the values of contact angle 4>q and are often the only experimentally obtainable quantities. Thus, the third component of energy may be written as The total energy then is, after nondimensionalizing and dropping all *'s. By scanning through the range of possible solutions the total energy E can be computed as a function of 0Â„ the contact angle which is contained in the solution for meniscus shape, including that corresponding to the specified value 4)^. The decision regarding stability of the meniscus can then be arrived at based on whether the specified value of 4>o minimizes E or not. Â£3 = Y/^ cos4)^ sir] (3.11) (3.12) + j27iy(lH//)'/2^ + 7ir;cos
PAGE 79

62 3.4 Numerics The Differential equation governing the meniscus shape, the LaplaceYoung equation, is nonlinear and has to be solved numerically. The integration of the nonlinear ordinary differential Eq. (3.4) is performed using the fourth-order Adams predictor-corrector method, using the Runge-Kutta method for starting. Onehundred and one points were used along the zdirection, as deemed sufficient from grid retlnement studies. The procedure adopted was to use the shooting method beginning from the fixed base at the edge of the die. Shooting was performed at various starting angles, henceforth 0b Â• The resulting slopes at the top of the meniscus were then stored. For any specified slope at the top, as required for a given contact angle, the profile was determined by using a bisection method by sweeping through the solutions stored from the initial value treatment. The integrals in the expression for energy were computed using the trapezoidal rule. 3.5 Results and Discussion As previously mentioned, the meniscus profiles are significantly influenced by Bond number (Bo), applied pressure (Ap), aspect ratio (h^ /T^,), and the value of contact angle. Our procedure was to specify the Bond number and aspect ratio, thereby obtaining the value of rb from Eq. (3.3) , and meniscus height from aspect ratio and rb. As already mentioned, an interesting feature of our solutions for this model of Edge-defined Film-fed growth is the presence of multiple solutions for a given

PAGE 80

63 contact angle. Unless a methodical sorting is performed, a numerical solution procedure may not pick up all the possible solutions. In the work of Katoh et al. (1990), multiple solutions seem to have been obtained only in certain cases. This may be due to their operating regime of parameters, such as Bond number, aspect ratio etc. From our solutions it appears that the tendency to obtain multiple solutions and interesting profile shapes is greater for taller menisci, and is also influenced by the applied pressure. 3.5.1 Basic mMeniscus Behavior with Upward Pulling Consider a typical case (Case 1) shown in Figure 3-2. The parameters are Bo = 4x10 "^ h,/rb (aspect ratio) = 0.5, Ap = 0. The curves corresponding to 0^ against 4)^,, and r,/rb against (t>^ clearly show the multiple solution behavior. Thus, there is only an interval of possible contact angles 0c that can be achieved at the top of the meniscus. Outside this range of angles no solutions exist. This is not due to instability of possible equilibrium profiles, but because it is not even possible to obtain profiles satisfying the LaplaceYoung equation. Also shown in Figure 3-2 are sample meniscus profiles, corresponding to (f>o of 60Â° and 90Â°, respectively. As shown in Figure 3-2, when the value of 0o is specified, and E is calculated, distinct extrema are obtained at specified values of o. Among the multiple solutions, the one corresponding to a minimum on the curve is the physically realizable stable meniscus profile. When either a maximum or a non-extremum point arises, we may surmise that such solutions belong to the unstable branch.

PAGE 81

64 That the Young's condition is in fact the condition for stability is straightforward to show, say in the case of a sessile drop. In such a (constant volume) case, when gravity is ignored, the drop assumes the shape of the arc of a circle (Dyson 1978). Then the minimization result becomes merely a geometric condition compatible with minimum surface area. When gravity is considered the drop no longer remains the element of an arc , since minimization of surface free energy is not sufficient, it being necessary to minimize the total energy including gravitational potential. In the case of menisci above infinite baths, Pitts (1976) and others have shown the minimization principle via variational methods. It is the conventionally held view that Young's contact condition can be obtained both from consideration of forces as well as energy (however, see Johnson (1959) and discussion therein for criticism of this view). However, it is interesting to see that in the present configuration, multiple shapes corresponding to energy minimum, maximum and even non-extremum may all be observed at the given contact angle 4>o. Furthermore, solutions to the LaplaceYoung equation, which is after all an equilibrium relation across the interface, are all equilibrium solutions. From the standpoint of local energy minimization, there may be more than one stable shape. Thus the end condition and the path of development of the meniscus selects the stable solution. However, in a process in which the aspect ratio is varied, an unstable meniscus may form first. This can happen if the new trijunction location corresponding to this shape is closer to the initial trijunction location. Then, in achieving the stable solution profile for the given height, one may pass through an unstable solution profile.

PAGE 82

65 As discussed in connection with the non-dimensional form of the LaplaceYoung equation, for low Bond numbers, the term corresponding to gravity is small compared to the curvature terms on the left hand side of Eq. (3.1). Thus the Bond number serves principally to relate the dimensions of the meniscus to the capillary length scale according to Eq. (3.3). However, for higher Bond numbers, the hydrostatic term on the right hand side of Eq. (3.1) becomes comparable to the curvature terms on the left hand side, and gravity does significantly impress on meniscus profiles, as will be illustrated in the results to follow. For taller menisci, i.e., higher aspect ratios, solutions to menisci appear to offer a wider variety of shapes and range of angles as well as a greater tendency to multiple solutions. Retaining Bo =4x10'^ and Ap=0, as in Case 1 , we can increase the height of the meniscus. For example, for hjr^, =2, which represents a taller meniscus, the available range of contact angles will shrink. Thus, as expected, when the height of the meniscus is increased, in the absence of pressurization, a large number of solutions to the LaplaceYoung equation fail to exist. Cases 2 and 3 (Figs. 3-3 and 3-4 respectively) are illustrative of the effects of pressure and Bond number on menisci with a higher aspect ratio where solutions with selected contact angles are shown. For Case 2, Bo=2.5xlO'*,Ap=5%, hjr^ =2. For this low Bond number case, from Figure 3-3, we notice that for a contact angle of 16Â°, three equilibrium profiles are obtained. The possible contact angles at the top were in the range of 18Â° to 26Â°. These profiles may not be seen in a physical situation, since the contact angles corresponding to the materials of the system may

PAGE 83

66 not lie in this range at all. Thus, ultimately the existence of menisci will be dictated by the range of contact angles ^ at the top permitted by the materials, the range of 0b for which the pinning condition is maintained at the edge of the die, and the stability of the equilibrium solution (of the LaplaceYoung equation) obtained. The contrasting situation to Case 2 emerges in Case 3 (Figure 3-4). Here the Bo = 2.5x10^ (h,/rb = 2,and Ap = 5% as in Case 2) is 10^ times that of Case 2. The available range of contact angles is considerably widened, and a greater variety of meniscus profiles than for the lower Bo in Case 2 can be obtained. Comparison of profiles with those in Figure 3-3 shows the difference that Bond number makes to the effect of pressure. Changes in curvature arise in the profile. This implies that for the given parameters, on the right hand side of Eq. (3.1), terms of pressure and gravity must compete to determine the shape of the meniscus. Thus, as mentioned before, for Ap 0, the Bond number can significantly influence meniscus shapes and permissible solutions. This can be seen from the scaling that we have employed in the non-dimensionalization. With fixed h, /rt, Eq. (3.3) shows that the height of the meniscus is scaled with respect to Bond number as h, Â« y/Bo. Hence, for given Ap, g and a, the meniscus height increases with Bond number. For higher Bo, by definition, the relative importance of gravity (Apgy) becomes more significant compared with the pressure term. Thus, the competition between gravity and pressure in the force balance leads to a variety of curvatures for higher Bond numbers. For lower Bo, a relatively small Ap can overwhelm the gravity term, and the curvature is less sensitive to location along the z-axis.

PAGE 84

67 Another aspect to note in Figs. 3-3 and 3-4 is that in the E vs. ^ curves both maxima and minima are obtained at specified Â„ = 130Â°, Figures 3-5 (c) and (d), the local minimum occurs at around 100Â° and there appears to be a nonextremum point corresponding to the specified indicating that although only one meniscus profile exists, it is not a stable one. Thus, for this case the decision regarding the stable meniscus profile becomes more involved. In summary, the numerical simulation reveals that depending on the range of parameters a variety of solutions are possible for a given contact angle. Multiple energy minima, implying multiple stable shapes, or no minima, implying unstable shapes may be obtained for a given contact angle. 3.5.2 Effect of Pull Direction on Meniscus Behavior The issues of the selection of a stable meniscus shape out of the possible solutions, and the implications of variation of h,/rb, as well as the direction of pulling on meniscus stability, can be better addressed through the energy considerations stated earlier. In studying the variation of meniscus shape and top radius with height the effect of the following factors is investigated. They are upward/downward pulling

PAGE 85

68 with respect to direction of gravity, Bond number (Bo), pressure difference between the liquid and gas phases (Ap), and aspect ratio (h,/rb). The results to be discussed in the following are obtained by fixing rb and varying hjr^ from 0.5 to 2.5 upward and downward, respectively. Solutions to the LaplaceYoung equation are obtained for each value of hjTf, and the range of available solutions in terms of (^<. and resulting radius at the top (expressed as r,/rb) are determined. Figure 3-6 compares the characteristics of the meniscus for the upward and downward pulling situations, with Bo = 10'^ and Ap = 0. The figures display the solution range of ^ versus h<./rb, the solution range of r^rb versus h,/rb, and possible meniscus profiles for a prescribed contact angle ^ = 10Â°, shown in the bottom of Figure 3-6.

PAGE 86

69 With Bo = 10"^ when hjr^ is less than 1, the hydrostatic pressure is small relative to the curvature terms on the right hand side of the LaplaceYoung equation. Hence the meniscus profiles are essentially determined from the geometrical constraints alone, meaning that within this regime the specification of r^, at the bottom and ^ at the top dictates the meniscus shape. Figure 3-6 clearly indicates that the effect of pulling direction is noticeable only when the aspect ratio is high enough. The critical value of h./rt for the hydrostatic pressure to exert a significant effect is Bond number dependent. It appears that the more appropriate indicator of the importance of hydrostatic pressure is (h^/rJ^Bo, i.e., an effective Bond number defined according to h, instead of rb. We next focus on the interaction of the direction of pulling, level of pressurization and variation of h,/rb. Figure 3-7 depicts the solution characteristics for the case of Bo = 10^ Ap = 1%, i.e., there is an excess pressure inside the meniscus. Regarding the overall range of permissible meniscus solutions, both in terms of 0^ and r,, the direction of pulling does not make any noticeable difference, seemingly indicating that the externally imposed pressurization dominates the hydrostatic pressure. However, depending on the value of h,/rb the detailed meniscus shapes can be significantly affected by the pulling direction. For example, with h,/rb less than or equal to 1, little difference is found between meniscus shapes formed by upward and downward pulling. For h,/rb of 2 or higher, however, the differences can be substantial. Specifically, for h,/rb = 2, with upward pulling, three possible meniscus shapes exist, while with downward pulling, only one profile can be obtained.

PAGE 87

70 When h./fb is raised to 2.5, for either direction of pulling only one meniscus profile is obtained, but the location of trijunction point for the downward pulling case differs from that of the upward pulling case. The energy profiles corresponding to the cases in Figure 3-7 are presented in Figure 3-8; for all the values of h, /rt considered here, both directions of pulling yield energy profiles that appear essentially identical. The energy curve for h^ /t^, = 2 shows that for the upward pulling cases, all three possible meniscus shapes contain virtually the same amount of energy despite their noticeable differences in trijunction location. Significantly, with hJVi, of 2 or higher, the energy content for both upward and downward pulling cases at <^o=70Â° correspond to either local maxima or nonextrema. Thus, for either direction of puUing, the meniscus shape can no longer be stable for such high aspect ratios. For larger h,/rb, the E value at the specified 0o is a local maximum of the energy profile. Hence, all such menisci are not stable to small disturbances, and tend to depart from the given meniscus shape and contact angle. Since here the static contact angle does not provide a minimum energy state, the meniscus, assuming it does not altogether collapse, will tend to seek such a state by changing shape and contact angle in a dynamic manner, 3.5.3 Ouasi-equilibrium Motion with Hysteresis Contact angle hysteresis is generally acknowledged to be a consequence of three factors: (1) surface inhomogeneity, (2) surface roughness, and (3) impurities on the surface (Miller 1985). In this third set of simulations we investigate the quasi-

PAGE 88

71 equilibrium path negotiated by the meniscus, by incorporating a simplified hysteresis model at the top trijunction, i.e., on our modeled flat plate. Specifically, the top plate is oscillated about a mean height h/ and the assumption is made that the time scale of these oscillations is such that a succession of equilibrium states is achieved. In other words, the LaplaceYoung equation is solved for each position h<.(t) of the oscillating plate. The contact line at the die is allowed to remain pinned at the edge, as long as the Gibbs' inequality is not violated. At the top, the boundary condition is modified by the hysteresis model adopted for the top trijunction, which is defined by the following generally observed criteria: 1. The contact line does not move, i.e.,re(t) is a constant when the contact angle varies between and 4>f^, where ^ and are respectively the advancing and receding contact angles (Dussan V. 1979). 2. When the contact angle <^a or <^r is reached, contact line motion ensues. To simplify the situation, we assume that the contact angle does not vary with contact line motion. i.e.,c=^ or ^ depending on the direction of motion. 3. It is imposed that an advancing contact line cannot recede without passing through the hysteresis range 0^ to 0r. Thus, if the oscillatory motion of the plate induces an impending recession of an advancing contact line, our model would perforce pin the radius r,(t) at such a location and require 0, to travel from 0^ to 0r before allowing recession. Similar conditions hold for the receding angle, 4>^. Schematic presentation of the hysteresis model is shown in Figure 3-9. An additional simplification resulting in the model is relevant in the presence of multiple

PAGE 89

72 solutions. When the contact line is advancing, for instance, we choose that solution for which T^{t) is closest in the direction of impending motion. So also for a receding contact line. Thus, irrespective of the stability of the equilibrium profiles, proximity of the contact line is made the criterion for choice of a solution from those available. Actually, in the situation where dynamic contact angles exist, energy criteria in terms of the existence of a minimum of the energy curve ( E versus ) is no longer applicable. In particular, the wetting component of energy E^^atd, namely TTj^{y,i-y,g) cannot be replaced by an expression such as 7rre^7,gCos<^e. for otherwise (7,i-7,g) would cease to be a material property. Thus, in passing through the successive equilibrium states, the stability of menisci cannot be deduced in terms of an energy minimization principle. From a microscopic viewpoint, if we consider that the hysteresis range j, to <^r merely represents the apparent contact angle behavior, and that in actuality, the microscopic contact angle is (f)Â„ (which is the Young's contact angle, not related in any straightforward way to 0^ or then we are justified in replacing (7,r7sg) by 7igCOS(^oHowever, on a microscopic scale, the phenomenon of hysteresis is occasioned by the presence of asperities (Oliver et al. 1977, Eick et al. 1975) and surface chemical inhomogeneities (Neumann 1972) or by several other effects which are still a matter of investigation. Thus, the wetted energy expression will need to be modified to 7,gCos
PAGE 90

73 It may be pointed out that the situations presented in this section and the two previous ones represent two limiting types of behavior. In sections 3.5. land 3.5.2, the boundary conditions and stability criterion, in terms of the minimization of energy are applicable to a situation where the time scale of motion of the top contact line is greater than that required for the meniscus to achieve equilibrium. Here, on the other hand, it is supposed that the time scale of contact line motion is much shorter than that required for reaching an equilibrium state, the latter being determined only by incorporation of the dynamics of the interface as well as the fluid contained in the meniscus. Traditionally, the following formula has been employed in models describing growth of a crystal via a pulling process. The variation of radius at the top of the meniscus with pull rate and time is given by (Thomas et al. 1986, Tatarchenko and Brenner 1980), = -^) tan(4) 4)^ (3-13) where (f>-o is the deviation in contact angle from the value that leads to fibre growth with non-uniform diameter. We attempt, in what follows, to evaluate the applicability of this expression in the light of hysteresis. In the event of 0 being a constant, it is obvious that the variation of the radius with time follows the same waveform but has an opposite phase with respect to the variation of meniscus height, hc(t). In this respect, the presence of hysteresis introduces a non-linearity in the above expression, in that the value of 0 in turn depends on r<.(t) and its rate of

PAGE 91

74 change. However, if a pinning condition arises, i.e.,re = constant, then the above formula indicates that, provided Up ^ dh/dt, <^ =
PAGE 92

75 given by h,(t), causes a corresponding sinusoidal motion of the top radius, r<.(t). The r, variation is here out of phase with h,(t) as represented by Eq. (3.13). It is noted that the radius variation in Figure 3-10 for the imposed amplitude AhÂ„ is quite wide, assuming values from 0.62rbto 0.83rb. Introduction of a hysteresis range gives rise, in Figure 3-11, to a square waveform for r<.(t). For the hysteresis range imposed, 0R = 85Â° to f^=95Â°, the pinning condition comes into effect when 4>^ lies in the hysteresis range. This can be seen from Figs. 3-11 (c) and (d), in which, while
PAGE 93

76 decreasing, while the angle remains constant. This is evidently due to the change in meniscus curvature from concave, as in Figure 3-10, to convex. For a higher meniscus, h/ = 1.5, for the same 4)^ and Ap shown in Figure 3-13, the radius variation is also not very substantial in magnitude and the waveform is not identical to h<.(t), indicating the non-linearity that is introduced by <^(t). When hysteresis is present, and for the case of fairly large hysteresis, and tall meniscus. Figure 3-14 shows an interesting situation arising due to the presence of multiple solutions. Note that in the previous cases, only one solution existed in each case, while here 2 or 3 solutions are obtained for certain heights. The choice of solutions is made as explained above. The marked excursions of radius evident in the figure arise when the outer solutions seen in Figure 314(a), are allowed to be attained by the imposed conditions 1^3. The hysteresis model in this case causes a noticeable change in the contact angle behavior (Figure 3-14 (d)) in comparison to the static model. There appears to be a Bond number dependence of the phase relationship between drjdt and dh,/dt, which is not reflected in the relation, Eq. (3.13). Thus, Eq. (3.13) relating the radius with height of the meniscus, seen in the light of a quasi-equilibrium path, does not appear to exhibit the appropriate features under contact angle hysteresis. In the absence of hysteresis, the Bond number dependence of the phase is not captured. When the LaplaceYoung equation yields multiple solutions, the choice of profiles is not straightforward, and the particular selection procedure adopted indicates an abrupt dynamic behavior which may not be justifiably represented by a quasi-equilibrium situation, and the full fluid-dynamic

PAGE 94

77 consideration is required to assess the implication of such behavior. We also reiterate that no account has been taken of the stability of the computed equilibrium profiles. 3.6 Conclusions 1. For the boundary conditions adopted here, which pertain some crystal growth processes, namely a fixed base radius, and Neumann boundary condition at the top, solutions to the LaplaceYoung equation were obtained using a fourth-order numerical integration scheme. Multiple solutions were obtained for a given contact angle at the top. 2. The energies of the menisci were computed and stability of solutions decided based on an energy minimization concept. While in many cases the energy versus 4>^ curves did indeed yield a minimum at a specified equilibrium contact angle 00, in other cases local maxima or non-extrema were obtained. Depending on the operating history of the pulling process, it may transpire that an unstable meniscus may form first, as for instance, when the trijunction point of an unstable shape at a new hjTi, is closer to that of the previous equilibrium position. 3. Unless the quantity (h./rJ^Bo is large enough, hydrostatic pressure does not considerably influence meniscus shape, and the effect of pulling direction is therefore not significant. The effect of external pressurization, Ap, also depends on the Bond number. For higher. Bond numbers, pressurization can markedly widen the range of

PAGE 95

78 available contact angles. Also, for sufficiently tall menisci, pressure (Ap) and gravity (Apgz) interact to cause menisci to assume a variety of curvatures. 4. For modest values of aspect ratio h,/rb, there always exist one or more stable meniscus shapes, and Young's contact condition is compatible with the requirement of energy minimization. For higher values of h^rb, although equilibrium solution satisfying LaplaceYoung equation do exist, they may not be stable solutions. 5. When the Young's contact condition does not yield a minimum in energy content, then the subsequent non-equilibrium dynamics of the meniscus needs to be simulated. Such treatments have been reported in connection with other situations (Young and Davis 1987, Haley and Miskis 1981) but have not been attempted for the pulling process. 6. Quasi-equilibrium studies of the dynamics of menisci have been made. A simplified hysteresis model was employed to relate the contact angle at the top to the contact line motion. The conditions enjoined by the model enforce pinning of the contact line when the contact angle lies between 0r and ^, the advancing and receding contact angles. It is found that the relationship in terms of waveform and phase of r,(t) with respect to h,(t) cannot be captured by the formula usually employed in crystal growth simulations.

PAGE 96

79 Figure 3-1 Schematic of meniscus corresponding to edge-defmed fibre growth process.

PAGE 97

1 Â« . 80 Figure 3-2 Characteristics of meniscus formation for Bo = 4x10"' . AP = 0, hjT^ =0.5, and upward puUingSO (a) variation of <^), with angle 0^ . (b) permissible range of tJt^ with respect to 0, . (c) profile shapes for different angles ^ . (d) profiles corresponding to 0,, = 60Â° . (e) E profile v.s.0, , for the specified value of 0Â„ = 60Â° . (0 profiles corresponding to Â„ = 90Â° . (g) E profile V.S.0, , for the specified value of 0, = 90Â° .

PAGE 98

81 U.J 0.45 0.4 0.35 j (d) 0.3 0.25 ( i) 0.2 0.15 / (ii) 0.1 0.05 \ n 1 -2 -IJ 1 -OJ 0 OJ 1 U 2 OJ -2 -IJ -1 -OJ 0 OJ 1 ij 2 I (j) minimiini 0 20 40 60 80 100 120 140 (i) minimum 0 20 40 60 80 100 120 140 Figure 3-2 -continued

PAGE 99

82 Figure 3-3 Characteristics of meniscus formation for Bo=2.5xl0* , AP=5%, K^h =2, and upward pulling. (a) profiles corresponding to = 16Â° . (b) E profile v.s.d), , for the specified value of qi)Â„ = 16Â° . (c) profiles corresponding to 0Â„ = 20Â° . (d) E profile v.s.0, , for the specified value of 0Â„ = 20Â° .

PAGE 100

83 Figure 3-4 Characteristics of meniscus formation for Bo=2.5xI0-^ , AP=5% K^h -2, and upward pulling. (a) profiles corresponding to 0^ = 60Â° . (b) E profile v.s.^ , for the spÂ°ecified value of 0 = 60Â° (c) profiles corresponding to 0^ = 90Â° . (d) E profile v.s.(^, , for the specified value of 0 = 90Â°

PAGE 101

84 Figure 3-5 Characteristics of meniscus formation for Bo=4xlO"' , AP=5%, hjTf, = 5, and upward pulling. (a) profiles corresponding to = 90Â° . (b) E profile v.s.0, , for the specified value of 0^ = 90Â° . (c) profiles corresponding to = 130Â° . (d) E profile v.s.0, , for the specified value of 0Â„ = 130Â° .

PAGE 102

85 (i) upward pulling (ii) downward pulling Figure 3-6 Effect of the direction of pulling on the meniscus characteristics with Bo=10"^ and AP = 0. (a) permissible range of contact angle. (b) permissible range of trijunction location. (c) permissible meniscus shapes at 0Â„ = 10Â° . (Note that with given Bond number" and hJr, less than 1 the meniscus shape is dominated by the geometrical constramt from the boundary conditions. .Hence the meniscus shapes are essentially the same between two pulling directions. The effect of pulling direction becomes noticeable only when h /r^ becomes larger.)

PAGE 103

86 Figure 3-7 Effect of pressurization on the meniscus shapes with Bo = 10"and AP = 1%. (excess interior pressure). (a) permissible range of contact angle. (b) permissible range of trijunction location. (c) permissible meniscus shapes at = 70Â° . (Note that the hydrostatic pressure effect becomes noticeable only for a higher h./r^ , as evidenced by difference caused by pulling directions.)

PAGE 104

87 (a) h^r, = 0.5 0.05 0.045 0.04 0.035 i 0.03 f E 0.025 \ 0.02 0.015 0.01 0.005 0 opwtfd puilins o o o aownwara pulling fii) noiv<]nremuin I (i) mininiism 30 40 50 60 70 80 90 100110 120130 c (c) hJT, = 2.0 40 45 50 53 60 65 70 75 (b) Vl, = 1.0 upwml pulliOE > downwara pulimg (ii) nofMxmmam (i) minimain I 20 30 40 50 60 70 80 90 (d) hJt, = 2.5 Figure 3-8 Energy profiles for Bo = 10"^ , AP = 1%, = 70Â° and various hJZf, . Both upward and downward pulling cases are shown.

PAGE 105

88

PAGE 106

89 0.60 0.40 1 0 40 75 115 150 time Meniscus behavior at lower Bond number, with no hvsteresis. r Bo = 10-^ AP = 5% h =0 5 rA,^^_^jooo^. ' (a) Meniscus profiles obtained for the duration of oscillation. (b) Variation of height h, with time. (c) Computed value of top radius r, versus time.

PAGE 107

90 c Fieure 3 11 Meniscus behavior at lower Bond number, with small hysteresis ( Bo = in-^ AP ^ h' _^J]JUiÂ„__SlZ^ = 85Â°) (a) Meniscus profiles obtained for the duration "^niation (b) Variation of height h, with time. (c) Computed value of top radius r. versus time (d) Computed variation of contact angle 0, versus time

PAGE 108

0.63 40 75 115 150 time Figure 3-12 Meniscus behavior at higher Bond number, with no hysteresis. ( Bo = 2.5x10-' AP = 5%. h ' = 0.5.6 . = = 100Â°V (a) Meniscus profiles obtained for the duration of oscillation. (b) Variation of height h, with time. (c) Computed value of top radius r, versus time.

PAGE 109

1.89 1.14 I 0 40 75 115 150 time Meniscus behavior at higher Bond number, with no hysteresis ( Bo = 2.5x10.AP = 5% h / = 1.5.6 . = 6. = innÂ°S (a) Memscus profiles obtained for the duration of oscillation (b) Variation of height h, with time. (c) Computed value of top Radius r, versus time.

PAGE 110

93 hc/rt, 1.14 40 75 115 150 time 92 (d) In/ i,y 1 / 0 40 75 115 150 time Figure 3-14 Memscus behavior at higher Bond number, with hvsteresis fÂ° = = J l.l^L5^,.^I10!^.^j00!}. a Memscus profiles obtained for the duration of oscillation (b) Variation of height h, with time. (c) Computed value of top radius r, versus time. (d) Computed variation of contact angle 0, versus time

PAGE 111

Figure 3-15 Meniscus behavior at higher Bond number, with small hysteresis range. ( Bo = 2.5x10-' AP = 5%. h,' = 0.5. 6 , = 90Â°. 6^ = 88Â°') . (a) Meniscus profiles obtained for the duration of oscillation. (b) Variation of height h, with time. (c) Computed value of top radius r, versus time. (d) Computed variation .of contact angle versus time.

PAGE 112

CHAPTER 4 NUMERICAL TECHNIQUES FOR SOLVING FREE AND MOVING BOUNDARY PROBLEMS WITH APPLICATION TO THE EFG PROCESS 4.1 Introduction As already discussed, prediction of the interface shape and movement in the EFG process is possible only if a complete dynamic model incorporating all the geometric and thermal aspects of the growth process is available. In order to facilitate such simulations, various computational elements must be developed, including appropriate time stepping and spatial discretization schemes, thermal field solutions in the irregularly shaped domain by the melt/gas meniscus and the melt/solid interface and method for tracking the free and moving boundaries (Tseng et al. 1992, Imaishi et al. 1991, VoUer et al. 1990, Floryan and Rasmussen 1989, Brown 1988, Salcudean and Abdullah 1988, Viskanta 1988, 1985 and 1983, Crank 1984, Ungar and Brown 1984, Epstein 1983). A literature survey of some existing methods for solving the Stefan problems is given in section 4.2. Then the issues arising from developing the front tracking techniques are investigated in section 4.3. We assess the applicability of the oft-used quasi-stationary approximation for both small Stefan number and large Stefan number cases. A new generic interface tracking procedure, based on a combined L^grangian/Eulerian framework, is designed; 95

PAGE 113

96 this procedure has been developed in conjunction with the previously mentioned time stepping scheme and the meniscus model to simulate the EFG process. After the presentation of the necessary numerical techniques, the EFG process physics is described in section 4.4. Fluid, thermal, and geometrical related non-dimensional parameters in thin crystal growth via the EFG technique are evaluated. Based on the order of magnitude estimation, a mathematical model is developed for the energy transport of the EFG process in section 4.5. Combining the front tracking and moving gridding techniques, time stepping scheme, solution procedure of the meniscus shape and contact condition, and the detailed thermal fields in both liquid and solid phases, a complete numerical tool for the EFG process is presented in section 4.6. The model is then applied to analyze cases with two Stefan numbers, i.e., St = 1 and 0.024. The steady state solutions corresponding to the specified boundary condition are obtained for both cases. It is identified that an appropriate scaling which accounts for the relative magnitude of the latent heat is necessary to properly capture the interface movement and to improve the computational efficiency. Based on the computed results of the two cases with different values of Stefan number, some conclusions for steady state solutions are made in section 4.7. The obtained steady-state solutions will be used as base solutions to study the dynamic behaviors of EFG system in Chapter 5.

PAGE 114

97 4.2 Literature Review The characteristic feature of any phase change problem is the coupling of thermal fields in the different phases with free and moving boundaries that not only separate each phase, but also dynamically evolve. The free surface and propagating phase front make such problems nonlinear. Because of the nonlinearity, only a few exact solutions have been developed, and these are limited to simple geometries and boundary and initial conditions. For most practical cases, the solutions are obtained using numerical methods. In what follows a brief survey of the existing methods for solving the phase change problems is presented. 4.2. 1 Exact Solutions One-dimensional heat transfer with phase change was first posed mathematically by Stefan and later solved analytically by Neumann (Crank 1984, Ozisik 1980, Carslaw and Jaeger 1959). The exact solutions are limited to a number of idealized situations involving semi-infmite regions or infmite regions subject to simple boundary and initial conditions. To date, no exact solutions have been developed for phase change problems where the phase change front is a function of multiple space coordinates. 4.2.2 Approximate Solutions The integral method, which dates back to von Karman and Pohlhausen, who used it for the approximate analysis of boundary layer equations, was applied by Goodman (1958) to solve a one-dimensional transient phase-change problems.

PAGE 115

98 This method provides a relatively straightforward and simple approach for approximate analysis of one-dimensional transient phase-change problems. 4.2.3 Numerical Solutions Many kinds of crystals grown by a EFG technique usually have high melting temperatures; for example, sapphire (AI2O3) with a melting temperature of 1046 K. Therefore, the need of simulating the solidification process has led to a development of several numerical methods. Most of them were summarized in Crank (1984), Wilson et al. (1978), and Ockendon and Hodgkins (1975). Chemouskou (1970) proposed the isotherm migration method, which interchanges the dependent variable with the space variable for one-dimensional phase change problems. Crank and Gupta (1975) were the first to apply the isotherm migration method to phase change problems in two dimensions. The advantages of this technique, in tracking a moving isotherm at the fusion temperature, are obvious, since the phase boundary is itself an isotherm, provided the phase change takes place at a constant temperature. Results predicted by this approach have been found to be satisfactory. Saitoh (1978) extended the twodimensional isotherm migration method from regular two-dimensional geometries to arbitrary shaped, doubly connected two-dimensional geometries. The moving heat source (or the integral equation) method, originally applied by Lightfoot (1929), is based on the concept of representing the liberation (or absorption) of latent heat by a moving plane heat source (or sink) located at

PAGE 116

99 the solid/liquid interface. The temperature at any point will be due to the superposed contributions from this moving source (or sink). In this method, the analysis of the phase-change problem is transformed to the solution of an integral equation for the location of the liquid/solid interface. Ozisik (1978) discussed the use of a moving heat source in the conduction equations to account for the latent heat effect. He presented a mathematical development that explicitly casts the moving boundary problem into a standard heat conduction problem with a moving heat source. By doing this, the solution can immediately be written in terms of Green's function, and numerically implemented. In general the numerical methods for solving the phase change problems can be divided into two categories: the fixed-grid approach and the transformedgrid approach. The fixed-grid approach implies that a grid is fixed in space and the interface conditions are accounted for by the definition of suitable source terms in the governing equations. In the transformed-grid approach, the governing equations and their boundary conditions are cast into a generalized curvilinear coordinate system, so the grid might adapt with the moving freezing front. A comparison between transformed grids and fixed grids on a test problem of the melting of gallium was made by Lacroix and VoUer (1990). They have found that if an efficient grid generation can be achieved, the CPU usage of the transformed grid is very close to that of the fixed grid. Furthermore, a fixed grid can produce accurate predictions with the same order of mesh size as that used by a transformed grid. Thus, according to these authors, no conclusive preference

PAGE 117

100 can be made for either method. However, the tests covered only a limited range of grid and time resolution; more effort in this direction will be needed to further clarify this issue. One well known method to deal with phase change problems on fixed grids is the enthalpy-porosity technique. This approach easily permits solutions for substances whose change of phase occurs over a range of temperature. Brent et al. (1988) have adopted this technique in conjunction with the SIMPLE procedure proposed by Patankar (1980) to simulate the two-dimensional convective melting of pure gallium in a rectangular cavity. The enthalpy-porosity technique was validated by comparing its results with the experimental results obtained by Gau and Viskanta (1986). A review of available implicit difference enthalpy-porosity schemes can be found in Voller (1985). One of the drawbacks of the enthalpy formulation, particularly for the case where the change of phase occurs at a single temperature, is that a plot of temperature against time for a given node tends to exhibit a "plateauing" tendency. This arises out of the fact that a computational grid models or represents a discrete region in space. Obviously, it requires a finite amount of time to melt a discrete region. As a consequence of a node being held fixed at a single temperature for a discrete amount of time, the effect is also felt by the surrounding nodes, causing a plateauing of the temperature. Voller and Cross (1981) have developed a smoothing technique that can be applied to a final set of numerical results. This smoothing technique has the effect of bringing the time-

PAGE 118

101 temperature history of a given node into better agreement with other published results. Treatments of the phase change and the associated transport process based on the fixed grid method have been developed to treat quite complicated problems. Bennon and Incropera (1987) have used a volume averaged continuum model with the SIMPLE algorithm (Patankar 1980) to investigate solidification of binary, aqueous ammonium chloride solution in a rectangular cavity. Beckermann and Viskanta (1988) have also used the same numerical method to solve the complex phase change problem. Shyyand Chen (1991,1990) have combined the enthalpy formulation and computational scheme based on an adaptive grid in curvilinear coordinates to solve the phase change problems with high Rayleigh and Marangoni numbers, and under both normal and reduced gravity conditions. With relatively finer and improved resolutions, various detailed transport phenomena and interface shapes have been reported. Recently, a turbulence model has been incorporated into this framework to predict the solidification aspects of a continuous ingot casting process (Shyy et al. 1992a), yielding a highly advanced tool for analyzing transport processes during solidification. Crowley (1983) applied an explicit finite difference scheme to an enthalpy formulation and coordinate transformation techniques to study the crystal pulling by the Czochralski technique. She applied her model to analyze the initially steady Czochralski system with a base temperature perturbed at the melt inlet.

PAGE 119

102 Based on the computed response of trijunction location, she concluded that the simulated configuration is unstable to such disturbances. A numerical solution for melting around a vertically heated cylinder incorporating the effects of natural convection was obtained by Sparrow et al. (1977). They used the coordinate transformation to immobilize the solid-liquid interface and, consequently, appear to have been the first investigators to apply it to a multidimensional convection-dominated Stefan problem. However, the higher order curvature terms for the interface were neglected, limiting the validity of the method to cases where the radius of the interface varies only slightly. Furthermore, a quasi-steady assumption was invoked where the effects of interface motion were ignored; hence the method can treat the problem only with very modest phase change rate. Thompson and Szekely (1989) have applied this quasi-steady approach to predict the phase change problem in a cavity. The above mentioned works are based on the finite volume/difference formulations. Silliman and Scriven (1980) and Cuvelier and Schulkes (1990) have used the finite element method (FEM) to calculate the predictions of steady viscous free surface flows. In their analysis, only liquid/gas interface is considered, and phase change is not included. Albert and O'Neill (1986) and Lynch (1982) used a finite element method with deforming elements to track a moving phase front and to model phase change problems. Ettouney and Brown (1983) have used a finite element analysis to model the steady EFG process. Several variations of the methods have been developed

PAGE 120

103 by this group (Sackinger et al. 1989, Derby and Brown 1986), depending on the choice of coordinates (original or transformed), distinguished boundary condition (local heat balance or thermal equilibrium), and the method of solving the resulting system of equations (successive approximation method or Newton's method). Their numerical method was limited to the steady solidification problems. Crochet et al. (1987), Dupont et al. (1987) and Wouters et al. (1987) investigated the oscillatory fluid flow in a horizontal Bridgman technique. The method was based on an implicit time-marching technique using FEM but was not applied to the moving solid-liquid interface problem. Other works in this area can be found in (Heinrich et al. 1989, and Dantzig 1989). In order to accurately resolve the solid/liquid interface for the isothermal phase change problem encountered here, the explicit front tracking method is adopted. A new technique has been developed, as presented below. 4.3 Front Tracking and Moving Grid The focus of this section is the numerical simulation of interface motion during melting/solidification of pure material. First, we assess the applicability of the aforementioned quasi-stationary approximation (Lacroix 1989, Thompson and Szekely 1989, Crowley 1983, Sparrow et al. 1977) for tracking the interface motion. Such an approximation for small St as well as large St is numerically examined. Next, a generic interface tracking procedure is designed, which

PAGE 121

104 overcomes restrictions of single-validness of the interface imposed by commonly used mapping methods (Lacroix 1989, Thompson and Szekely 1989, Crowley 1983, Sparrow et al. 1977). This method incorporates with ease interface phenomena involving curvature, which assumes importance at the smaller scales of a deformed interface. Among the methods proposed in the literature to track moving fronts (Floryan and Rasmussen 1989, Crank 1984) one commonly used approach (Lacroix 1989, Thompson and Szekely 1989, Crowley 1983, Sparrow et al. 1977) involves coordinate transformation to map the irregularly shaped physical domain onto a regular computational domain. The solutions may then be obtained in the transformed domain and the interface position updated. Often a further simplification is made by assuming that the time scale for interfacial motion is long compared to thermal relaxation times. This justifies dropping terms representing grid movement from the equations governing interface motion. This procedure is termed quasistationary and is appropriate for slow-moving interfaces. 4.3.1 Assessment of the Quasi-stationary Approximation As a result of the nonlinearity of the equation introduced by the motion of the phase change boundary there are but a few exact solutions of problems concerned with melting and solidification (Crank 1984, Ozisik 1980, Carslaw and Jaeger 1959). One such problem is the melting of a solid (Fig. 4-1) confined to a

PAGE 122

105 semi-infinite region. The solution is known as Neumann's solution (Crank 1984. Ozisik 1980, Carslaw and Jaeger 1959). We will use it to verify the developed numerical algorithm. For a pure conduction problem with phase change between liquid and solid, the governing equations, in dimensional form, for the energy transfer and interface movement are respectively (Crank 1984) The partials with respect to n represent derivatives in the direction of the local normal to the interface. The following non-dimensionalizing procedure is adopted. X=x/lr , Y=y/lr , S=s/lr , N=n/lr , 1^ being a representative length scale to be chosen for the specific problem studied. r'= ol-^\^, p* = Pj /pi and 6, = (Tj Tm)/(T, T^ , where Tn, and T, are, respectively, the melting temperature and the imposed temperature at the appropriate boundary. Consider heat transport in the liquid phase only, i.e.,^, = 0. Then, Eqs. (4.1) and (4.2), discarding the superscript, become: i = l,s (4.1) (4.2)

PAGE 123

106 dd _ d^e . d\ (4j) + dr' dX^ BY^ _ dd _ p * dS (4.4) 'dN IFTr where St = Cp( TÂ„-TJ /If is the Stefan number. If X=S(Y,t), then Eq. (4.4) specialized to an isothermal interface can be rewritten as Â— = Â£lÂ— ^ (4.5) dY dX St dr' Following the transformation from Cartesian(X,Y,r) to curvilinear (^,v,t') coordinates (Thompson et al. 1985), such that X = X (|,r/,T*),Y = Y{^,7],t) and r* = T, Eqs. (4.3) and (4.4) become: d + 1(X Y -YX )d, + -(Y^ -XJ )d (4.6 a) ^[(Ye\ (Y^d)^ [1 ^ (1[-(X5), ^ iX^S)^y ] = -A[^+I(xy YX)s, + l(y^ A:,y)5j St dr J " ' " y ^ ' ^ ' (4.6 b) where (4.7)

PAGE 124

107 Eqs. (4.6 a) and (4.6 b) constitute the complete set of the equations in the generalized coordinates for energy transfer and interface movement. In the course of the computations, at every time step, both Eq.(4.6 a) and (4.6 b) need to be iteratively updated until they are simultaneously satisfied everywhere in the domain. Standard second-order central differences are used for all spatial derivatives. The quasi-stationary approach adopted by many researchers, e.g., Lacroix (1989), Thompson and Szekely (1989), Crowley (1983), Sparrow et al. (1977), simplifies the above equations by dropping all the terms involving coordinate movement, i.e.,XÂ„ Y^. To test whether such a simplification is acceptable, we design a situation where the interface moves as a planar front. In addition, the temperature in the solid is considered uniform, i.e., 6^ = 0. Thus the heat transport is effectively one-dimensional and occurs in the liquid phase only. Boundary and initial conditions are therefore given by Crank (1984), B.C. e 1, X=0, T>0 (4.8) e 0, X>5, r>0 I.e. erfi X (4.9) 5(X,rÂ„=0.1) = 2\f^,

PAGE 125

108 where X is the root of the following equation Xe'^'em = A (4.10) The exact solution (Crank 1984) for Eqs.(4.3) and (4.4), with p* = 1, with the boundary and initial conditions given in Eqs.(4.8) and (4.9) is '^77=^ (4.11) em Two numerical solution procedures are designed to compare the solutions obtained from the complete form of Eqs. (4.6 a) and (4.6 b) and the quasistationary form. These are, (1) Full treatment: A standard procedure involving backward Euler time stepping along with a second-order central difference spatial discretization scheme is used to solve Eqs. (4.6 a) and (4.6 b). This full set of equations is solved iteratively in a coupled manner to continually update the nonlinear coefficients resulting from the coordinate movement and transformation. (2) Quasi-stationary treatment: By invoking the quasi-stationary assumption, coordinate movement terms are neglected. On dropping the grid movement terms, the equations governing energy transfer and interface movement are decoupled and only a simple explicit procedure is needed to update the interface location. Three values of the Stefan number were chosen to test the performance of the above two numerical approaches. Eleven grid points are used along the x-

PAGE 126

109 direction in each case. For small St. i.e.,0.1303,as seen in Fig. 4-2 (a), both full and quasi-stationary solutions are in close agreement with the exact solution in terms of the interface trajectory as well as the temperature profile, although the full approach is marginally superior. However, as expected, with increasing St, the quality of solutions obtained from the quasi-stationary method is progressively degraded. In contrast, good agreement is maintained between the solutions of the full approach and the exact solution. Obviously, with a low St, the interface velocity is modest, and neglecting grid movement terms does not impact significantiy on the numerical accuracy. As St becomes larger, however, the simplification to the quasi-stationary approach is no longer acceptable. 4.3.2 A General Procedure for Interface Tracicing As mentioned in previous subsection, Eqs. (4.6 a) and (4.6 b) apply only to an isothermal interface. Also, the construction of those equations precludes the possibility of capturing branched interfaces. As a first step towards overcoming this difficulty, an alternate interface treatment is designed as follows. Consider the equation for interface advance k dd s s k, BN (4.12)

PAGE 127

110 The local normal to the interface is given by (4.13) N = |VF| dx where F = F(X,Y,t') is the curve defining the interface. The interface shape is defined in a piecewise fashion to facilitate handling of branched interfaces. Here a quadratic polynomial fit is performed for three successive nodal points at each point of the interface. Thus, at the i* point on the interface we designate the curve, Yi = a^X;^ + b^X^ + q , i.e.,Fi = Yi ( a^Xj' + bjX; + C; ) defines the interface. The a^, bj and Cj are determined from the known values of (Xj,Yj), j = i-l,i , i+1. The choice of the piecewise approximating polynomial is not restricted to the parabolic form. In fact, using circular arc elements is more convenient for handling multiple-valued interfaces and is in use in ongoing work. The local curve definition yields the derivatives (F,, Fy and F^ ) at each point on the interface. We may write Eq.(4.12) as In computing the interface normal velocity then, one seeks to obtain the derivatives F^, Fy and 6^ , dy. The derivatives of temperature may be obtained in the transformed coordinates itself. F^ and Fy of course are directly available in the physical domain from the curve fit. Thus the new coordinates of the interfacial points are obtained from: |VF| [ (FA F^^), (FA F^,). ] (4.14)

PAGE 128

Ill Xn^i ^X" ^ Â—-!^8t (4.15 a) dX |VF| y*' = y+ K-!^8t (4.15 b) BY |VF| where 5r is the non-dimensional time step size. Having obtained these new coordinates of the curve, the thermal field is solved for once again, the curve fit is h .... performed at the interface and a fresh interface position is obtained from Eqs.(4.15 a) and (4.15 b). All these procedures are performed in a fully coupled manner involving interaction among the temperature field, interface motion and grid movement at each iteration. This alternative interface ti^eatment is compared first with the results of section 4.2.1, where the thermal field is active in the liquid region only. The boundary and initial conditions are given by Eqs. (4.8) and (4.9). From Fig. 4-3, it is evident that for all three Stefan numbers tested, namely, 0.1303, 1.2216 and 2.8576, the accuracy of the present scheme is satisfactory and comparable to the full treatment in section 4.2.1, Fig. 4-2. As can be seen, the computed and exact temperature fields are in excellent agreement. The metiiod has been extended to the case where temperature gradients exist in both liquid and solid, and applied to study the development of a morphologically unstable interface (Shyy et al. 1992 b). 4.3.3 Assessment of Time-stepping Schemes for Phase Change Problems As already mentioned, the EFG process operates in the conduction dominated regime. Hence, the interaction between the time stepping and the

PAGE 129

112 convection schemes discussed in Chapter 2, while generically important, is not critical for our main interest here. In what follows, a comparative study is conducted to assess the relative performance between the backward Euler and Crank-Nicolson schemes for solving one-dimensional Stefan problems. The governing equations of the unsteady heat conduction problem in one phase, along with the energy balance requirement at the interface boundary, boundary conditions and initial condition are given by Eqs. (4.3), (4.4), (4.8) and (4.9). The computations start with the exact solution at t = 0. 1 as the initial condition, and stop at r = 0.2 with four different values of At with St = 0.1303. The predicted interface locations at r = 0.2 along with required CPU time using Micro Vax 3100 are listed in Tables 4-1 (a) and 4-1 (b). Two different spatial grid sizes. Ax = 0.0158 with 11 grid points and Ax = 0.0079 with 21 grid points, along with several time step sizes have been used in the assessment. Based on the Table 4-1 (a) the accuracy of the Crank-Nicolson method is very comparable to the backward Euler method. In terms of the CPU time requirement, both methods are essentially the same, as expected. It should be pointed out that in the interface updating procedure via the integration of the paths of marker particles, discussed in section 4. 3. 2, a first-order one-sided scheme is used to estimate the thermal gradients at the interface. Consequently, the marker velocity is first-order accurate, and the resulting interface location is not dictated by the choice of time stepping scheme. It is apparently this practice that

PAGE 130

113 is responsible for the comparable performance between backward Euler and Crank-Nicolson, However, as will become clear, due to the nonlinearity of the phase change problem, only a relatively small time step size, around 10'\ can be used to sustain a stable computation. Accordingly, there is little preference between the two choices of time stepping schemes. For all the computations to be presented in the following, the backward Euler scheme has been adopted. 4.4 EFG Process Physics In the following, the basic process physics will be first introduced along with assessments of the key control parameters including the Fourier, Bond, Rayleigh, Marangoni, Peclet numbers and aspect ratio. Then the appropriate governing equations and boundary conditions will be presented, and the numerical aspects needed for simulating such a process discussed in the next section. The basic EFG operation considered here is shown in Figures 1-1 and 1-2, including a schematic illustration and a description of the boundary conditions and notation adopted. A major feature of the present process is the relatively small size of the meniscus and the fibre diameter. As indicated by the Bond number, defined as Bo = (4.16) where rb is the lower contact radius on the die, 7,^ is the surface tension between the melt and the air, Ap is the density difference between the melt and the air.

PAGE 131

114 and g is the gravitational acceleration. With the size and properties given in Table 4-2, Bo is of the order of 10"*, meaning that the capillary effect is important in determining the meniscus shape. Three types of convection exist in the melt in such a system. The first one is the buoyancy driven convection caused by the interaction between density gradient and gravitational acceleration in the bulk of the melt. The second one is the thermocapillary convection caused by the surface tension gradient due to the temperature gradient along the melt/gas interface. The third one is forced convection caused by the flow of the melt towards the solidification front, due to the crystal pulling, in order to conserve mass. Based on the nominal values fibre diameter of 0.005 ~ 0.025 cm, and a temperature variation of 5 ~ 25 K from the die tip to the solid-liquid interface, the Rayleigh number is 0(10'*) for sapphire, according to the property values listed in Table 4-2. Since the interface separating the melt from the outside gaseous environment involves a large temperature gradient along the free boundary, the Marangoni effect should also be estimated. Again, with the reference values given in Table 4-2 and a surface tension of 0.06 dyne/cm-K for dy^/dT, the Marangoni number Ra = (4.17) Ma = (4.18) is about 0(1). The Peclet number measures the importance of heat transfer by

PAGE 132

115 convection, based on the pull speed, Up, in the growth direction relative to conduction Pe = ^ (4.19) is 0(10"'). For the present system, the small values of Ra, Ma or Pe indicate that the thermal convection within the melt is of a minor effect. Hence, neither thermocapillary nor buoyancy induced convection is included in the model; forced convection is included here mainly to account for the mass conservation constraint. Accordingly, in the present model, besides the transient and twodimensional conductional effects within the melt and solid, the only convection effect considered is the forced flow fluid caused by the continuous pulling of the crystal at the speed Up. 4.5 Mathematical Model In the following, a theoretical model of heat transfer in the EFG process under transient conditions is formulated. The governing equations of the heat transfer process in the melt and solid are presented along with the constraint controlling the shape of the melt/solid interface. The interface condition adopted here is based on phase change equilibrium (isothermal phase change condition) and local heat balance without considering the microscopic effects such as the Gibbs-Thompson effect or the anisotropic surface energy (Pelce 1988).

PAGE 133

116 We are interested in the axisymmetric situation, i.e. r-z coordinates, with which the crystal is grown from the melt. A schematic of the melt region of the present EFG model along with the associated boundary conditions is given in Figure 4-4. Table 4-2 gives the physical dimensions and the material properties of sapphire, while Table 4-3 summarizes the important parameters and their numerical values encountered in the present model based on Table 4-2. Both large St and small St cases are numerically investigated. The values of Bond number, Rayleigh number, Marangoni number, and Peclet number are all small here, and hence conduction is the dominant heat transfer mechanism within the melt and the crystal, and capillarity dominates static pressure in terms of controlling the shape of the meniscus. Under such a circumstance, the meniscus profile is close to a straight line (Kalejs et al. 1983a). Regarding the governing equations and boundary conditions relevant to the present system, the following presents the LaplaceYoung equation (controlling the shape of the melt/gas interface), energy transport equation, mass continuity equation, and movement of the solid/melt interface in the r-z axisymmetric coordinate system. (i) meniscus shape The LaplaceYoung equation governing the equilibrium meniscus shape in the axisymmetric geometry is

PAGE 134

117 where 7,^ is the surface tension between the melt and the air, f describes the meniscus shape, and Ap is the difference in applied pressure between the liquid phase and the gas phase. The bottom edge of the meniscus is defined by the die radius, r^, , and the top edge is defined by the crystal radius, r, . (ii) energy equation Melt The energy equation in the liquid phase in dimensional form is dt bz r dr dr dz dz where the subscript 1 denotes the liquidus phase condition, is the convection speed along the axial direction and f and h define the meniscus shape in r and z direction, respectively. The boundary conditions are given as follows : (a) central line .^=0, r=0, Oizi/i(0) ("^-22) dr (b) the melt inlet r, = r,. 0^r^r,,z=0 (4-23) where Tb is the base temperature at the melt inlet, and rb is the radius of the contact line between the meniscus and die, and (c) along the side surface of the melt, heat is transferred to the surroundings by the combined conduction, convection and radiation

PAGE 135

118 dT -k' = hfT^ -r ) + eo[{T^ r -r % r=M, 0^z
PAGE 136

119 dT i=0, r=0, /j(0)izi/ ("^-28) dr (iii) mass conservation Because of the low Peclet number, convection does not affect heat transfer in a significant way. The velocity in the melt is determined using the onedimensional continuity equation: where v^is the advection speed in the melt, which is determined by the local radius. (iv) conditions at solid/melt interface Finally, at the melt/sold interface z = h(r,t), both equilibrium solidification and conservation of energy are satisfied T = J = T (4.30 a) Ism kf-Il pMu^\ (4-30 b) ' dn ' dn ' ^ ' dt " The various issues of the specification of an appropriate contact condition at the trijunction point have been discussed in Chapter 3. Conventionally, a constant contact angle, in the spirit of the Young's equilibrium contact condition, has been used (Sackinger et al. 1989, Derby and Brown 1986, Ettouney and Brown 1983, Boley 1975), namely, -S = (u -'-^)tan(0(r) <^^ (4-31) where r, is the radius of the trijunction point, h, is the height of the trijunction

PAGE 137

120 point, (^(t) is the instantaneous contact angle of the meniscus at the trijunction point, (i>o is the equilibrium contact angle. The above equations are then written in dimensionless form, by using the following reference length, time, and temperature, respectively [Length]: 1^ = h (radius of the die), [Time]: t, = h^/a^, (diffusion time scale), [Temperature]: AT = T^ TÂ„ Using the upper case letters to designate the nondimensional parameters in space, T* = t /t^ to designate nondimensional time, and T; = / a, , i = i, s p = p, /p, , h = h, / hÂ„ k = k, / k,, Up =\i,tJ a, V, = V, rb / a, = (T, TJ / AT, 9, = (T, TJ / AT, G,* = G, / (p,Cp,ATa, / rj The dimensionless energy equation and corresponding boundary conditions for the liquid phase become

PAGE 138

121 = lÂ±iR^)^Â±i^), 0
PAGE 139

122 e = 0 (4.40 a) where St is the Stefan number. The interface shape is determined by thermal equilibrium and local heat balance. The exact location of the trij unction point which defines the radius of grown crystal is decided by enforcing a constant contact angle at the upper trijunction point, which has been employed by several authors (Sackinger et al. 1989, Derby and Brown 1986, Ettouney and Brown 1983, Boley 1975), i.e., ^ = (t; -Â—1) tan(0(r) -<^^ (4.41) where Up is the pulling speed, <^(t) is the instantaneous contact angle of the meniscus at the trijunction point, and <^o the contact angle of the meniscus at steady-state. At the lower trijunction point, on the other hand, the contact radius is fixed by the sharp edge of the die tip. Rewriting the heat equation in general form for both liquid and solid phases, ^i(r,z,t), i = l,s: (Â— +v;.Â— ) = r.[-L-r_(/?_i) +_Z_(^)] (4.42) ar* *az ' /?a/? m dz dz where V, = V, and V, = Up . Employing the transformation from the (R,Z,r*) coordinates to curvilinear the (^,17,7) coordinates,

PAGE 140

123 R = R(^rj,T) , e = |(R,Z,r*), Z = Z(?,ij,t) , 7j = rj (R,Z,r'), (4.43) T T = T the above general equation after dropping the subscript "i" can be rewritten as where q, = Z,^ + R,^ q^ = Z^Z, + R^RÂ„ q, = + Rj', J = Z^R, Z,R^ , WÂ„ = V,R, , and W, = V.R^ Regrouping the above equation, and treating the grid movement induced convection terms and cross derivative terms as pseudo-source terms, we move them to the right hand side of the equation: where the spatial derivative terms are treated based on the second-order central difference scheme in the context of second-order finite volume formulation, and the unsteady term is by the backward Euler method. As already discussed, the computational strategy designed here utilizes an explicit front tracking technique, based on a combined Lagrangian and Eulerian approach, to locate and advance the moving boundaries between the crystal and the melt, as well as the new meniscus shape. In this procedure, the markers along (4.44) (4.45)

PAGE 141

124 the interface are used to define and advance the location of the interface in time; the movement is tracked via a Lagrangian framework. The velocity of each marker is determined by the balance between the difference of thermal gradients in the solid and in the melt, and the latent heat absorption as given in Eq. (4.40 b). Based on the newly defined boundary shapes, a grid will be generated again and the rate of grid movement evaluated accordingly. The energy equation in each phase will be computed in an Eulerian framework. All these procedures are performed in a fully coupled manner involving interactions among the temperature field, meniscus and interface motion, and grid movement at each iteration. Within each time step, the iterations continue until all the governing equations, boundary conditions, and interface constraints are satisfied within a specified convergence criterion; here, it requires that the summation over the domain of the absolute imbalance of the nondimensional fluxes in each computational cell to be smaller 5xlO"^ 4.6 Numerical Simulation of the EFG Process The grid distribution, generated by an algebraic procedure, transfinite interpolation technique (Thompson et al. 1985), is clustered in the regions close to the central symmetric line and the outside boundary of the crystal and melt. The control surfaces are determined with a power-law clustering scheme over half of the computational domain

PAGE 142

125 = ^^(4r^'"^ i =0,l,...,iVr (4.46 a) 2 N r 2. = !!^{Lr, j=o,u...,Nz (4.46b) where 2Nr + 1 and 2N^ + 1 are the total number of control volumes along Rand Zdirection, respectively. The exponents m and n control the rate of clustering in the Rand Zdirection, respectively. A value of 1.55 for m and 1 for n is used in this study, i.e. the clustering is assigned along the Rdirection only. The grid distribution is symmetric about the lines R = F(R,t)/2 and Z = H(R,r)/2. The grid with finer resolution in the region close to the outside boundary of the crystal and melt is found important for yielding good solution accuracy. The complete determination of the meniscus profile involves the solution of the LaplaceYoung equation constrained by hydrostatic and externally enforced pressure difference, Ap, the Bond number, defined in Eq. (4.20), the meniscus contact angle at the trij unction point, the meniscus height of the trijunction point HÂ„ and the die radius, R^. There are also issues of static versus dynamic contact angles, as examined in Chapter 3. One set of solufions obtained using the linear meniscus profile, i.e., in non-dimensional terms, R^ = 1, H, = 0.188,R, = 0.763, is chosen to compare the solution of LaplaceYoung equation of Bo from 2.8x10-^0 2.8xlO-\with Ap/Apgr, = 5% and 0%. The computed meniscus profiles are quite close to linear, as shown in Figure 4-5. With the

PAGE 143

126 straight line assumption, the meniscus volume does not differ by more than 2 % from that calculated based on the LaplaceYoung equation. This observation is consistent with the assumption made in Kalejs et al. (1983a). For these parameters relevant to the sapphire fibre production, and under the condition that the static contact angle controls the trijunction point movement, the use of a straight meniscus profile is a reasonable approximation. Both large St (Case 1: St = 1) and small St (Case 2: St = 0.024) are computed. Thermophysical properties and physical dimensions of sapphire (AI2O3), processing parameters and the corresponding non-dimensional operating parameters are given in Tables 4-2 and 4-3. As can be seen, for simplicity, constant properties are used for both solid and liquid phases. This aspect is not expected to affect the solution in a major way. Due to the symmetry of the problem, only half of the physical domain is modelled. For the entire physical domain encompassing both the crystal and the melt, effectively a total of 41x41 nonuniform grid is employed. Due to the complex nature of the problem, such a grid size already requires substantial computing resources to support the numerical simulations. In the present work, both engineering workstations, including IRIS Indigo, DEC 5000 and 3100, and the Cray Y-MP supercomputer have been used. For those workstations, it typically requires ten hours or more of CPU time to complete a simulation. The initial condition of both St = 1 and St = 0.024 is shown in Figure 4-6, and the corresponding steady-state solution is shown in Figure 4-7 and 4-8,

PAGE 144

127 respectively. Density and thermal diffusivity are taken to be the same for both the melt and solid phases. Figure 4-7 shows the isotherms of the steady-state solutions of St = 1 in both phases, along with the grid lines used to conduct the computation. The interface is flat in the inner region, and becomes more curved toward the outer region due to the heat loss to the ambient environment. It is noted that because the outer boundary of the solid and melt are of different slopes, the local curvatures of the isotherms in these two phases are different. The grid system and the steady-state isotherm distribution of St = 0.024 are shown in Figure 4-8; due to the small value of Stefan number, the thermal gradients in the melt and the crystal regions are of different magnitudes. Furthermore, it is interesting to note that with St = 0.024 the solid/melt interface is convex toward the melt, while in previous section, it is found that with St = 1 the interface is convex toward the solid region. In terms of the computational practices, we discuss two scaling procedures utilized to non-dimensionalize and, more critically, normalize the various terms present in the governing equations and their implications on the computational efficiency. As it turns out, care should be taken in choosing the appropriate reference quantities because for many phase change problems this choice has a major impact on the amount of computing cost required to conduct a simulation. It is noted that in the present problem, the two energy transport equations are applied in the regions separated by a moving interface, causing a nonlinear coupling between the two equations. Hence, in the course of solving these energy

PAGE 145

128 equations with the backward Euler time stepping scheme, while a von Neumann type of stability analysis indicates that the computation is unconditionally stable for each of the energy equations, the nonlinear coupling caused by the interface movement makes it only conditionally stable. Numerical experiments indicate that, as expected, the range of the time step size acceptable for a stable computation depends on the distribution of the grid points in each phase. Table 4-4 presents two different scaling procedures and the resulting nondimensional equations of energy transport and interface movement. As can be seen, the only difference between the two procedures is the velocity scale; in choice 1, a standard characteristic diffiisional velocity is used, while in choice 2, the Stefan number is included in addition to the diffusional velocity scale. The main motivation of the second choice stems from the observation that as can be easily deduced from Eq. (4.40 b), with choice 1, the non-dimensional speed of interface movement is scaled with St. Accordingly, in non-dimensional terms, for the low St cases, the interface moves at a correspondingly slow rate. As already mentioned, even with the use of an implicit procedure, the computation is not unconditionally stable. Numerical experiments have indicated that, in terms of the non-dimensional values, the time step sizes of the set of governing equations resulting from both choice 1 and 2 have very comparable stability restrictions. Hence, considering the fact that the time scale of choice 2 is larger than that of choice 1 by a factor of 1/St, the relative computing efficiency of the two scaling procedures depends highly on the value of the Stefan number. As an illustration.

PAGE 146

129 computations have been conducted using both scaling procedures and starting with the identical initial conditions for the low St case studied here. In both computations, the non-dimensional time step size is At = 10'^ , which is about 10 times of the nondimensional diffusional time scale based on the smallest grid spacing, namely (Ax)^Va , where (Ax)^ is the smallest grid spacing and a is the thermal diffusivity. For both scaling procedures, numerical experiments indicate that this value of At is very close to the stability limit of the implicit Euler time stepping scheme adopted in the present work. Figure 4-9 compares the paths of the computing histories with both scaling procedures from an identical initial condition toward the steady-state solution subject to the same boundary conditions. The Stefan number chosen in the present work is 0.024. Substantially different convergent behaviors have been observed between the two scaling procedures. As already mentioned, the differences of the convergence characteristic is largely caused by the different orders of magnitude of the non-dimensional interface velocity yielded by the two scaling procedures. With choice 1, the interface speed is of the order of St, while the choice 2, it is of the order of 1. Consequently, choice 2 needs a far fewer number of non-dimensional time steps to reach the steady state solution. With choice 1, however, the system appears to go through a long transient period before it can reach the steady-state. As a check, the converged steady-state solution obtained from choice 2 is substituted into the computation based on choice l,and this solution is immediately accepted as the converged solution,

PAGE 147

130 confirming that with both scaling procedures, the computations eventually lead to the same steady-state solution. However, the computing costs of the two scaling procedures differs by an order of (1/St). Obviously, choice 2 is appropriate for the low St computations. 4.7 Conclusions In this Chapter, a thermocapillary model has been developed to simulate the characteristics of the EFG process. The model is axisym metric, solving the energy equation in both solid and melt phases separately, and tracks the movement and evolutions of their interface explicitly using a combined Lagrangian/Eulerian method. The meniscus behavior is governed by the LaplaceYoung equation subject to a specified contact angle at the trijunction point. It has been established that for the thin fibre process, the value of various control parameters such as the Stefan number, Rayleigh number, Marangoni number, and Bond number are modest, and hence conduction is the primary heat transfer mechanism within the melt and the crystal, and capillarity dominates static pressure in terms of controlling the shape of the meniscus. The quasi-stationary approximation for interface is found to be inaccurate for capturing the phase boundary except for small Stefan numbers. With a full treatment, the temperature field, interface calculation and grid movement all interact over each interaction and time step. Such a computational procedure is applicable to an isothermal interface. A generalized interface motion technique is

PAGE 148

131 developed, which is not limited to a single-value or isothermal interface. Computations show that this method is both accurate and versatile. Comparison of the two scaling procedures for the interface movement shows that the characteristic velocity scale should be defined to ensure that the solid/liquid interface advances at a nondimensional speed of order 1 . Since the interface movement is primarily responsible for the nonlinearity of the phase change problem, a conventional choice based on the diffusion scale, choice 1, will yield a set of equations that requires the computing cost to increase by a factor of 1/St. It is also interesting to note that with St = 1 , the solid/liquid interface is convex toward the solid region, while with St = 0.024, the interface is convex toward the melt.

PAGE 149

132 Table 4-1 Assessment of time stepping schemes for phase change problem. (a) Comparisons of the predicted interface location at r = 0.2 using backward Euler and Crank-Nicolson time stepping schemes for Eqs. (4.3) and (4.4) with 11 and 21 grids. 11 grid (Ax = 0.0158) 21 grid (Ax = 0.0079) At Backward Euler Crank-Nicolson Backward Euler Crank-Nicolson 10-^ 0.223 0.222 0.227 0.223 10-^ 0.226 0.225 0.226 0.226 10^ 0.226 0.226 0.226 0.226 10-^ 0.226 0.226 0.226 0.226 (b) Comparisons of the required CPU time on Micro Vax 3100 using backward Euler and Crank-Nicolson time stepping schemes for Eqs. (4.3) and (4.4) with 11 and 21 grids. 11 grid (Ax =0.0158) 21 grid (Ax = 0.0079) At Backward Euler Crank-Nicolson Backward Euler Crank-Nicolson 10-^ 0:12.28 0:13.07 0:16.77 0:21.00 10-^ 0:39.84 1:12.24 1:09.79 1:13.68 10^ 4:54.87 4:5967 9:15.15 9:14.34 10-^ 49:49.27 49:48.20 1:31:27.93 1:31:15.27

PAGE 150

133 Table 4-2 Thermophysical properties and physical dimensions of sapphire (AI2O3) used in the present EFG simulation (Case 1: St = 1 and Case 2: St = 0.024). The material property of sapphire (AL O^') Case 1 Case 2 (St = 1) (St = 0.024) k, (W/cm-K) 0.1 0.1 0.1 p, (g/cm') 3.05 3.05 3.05 Cp, (J/g-K) 1.26 1.26 1.26 e, 0.9 0.9 0.9 k, (W/cm-K) 0.1 0.1 0.1 p, (gW) 4.00 3.05 3.05 Cp, (J/g-K) 1.26 1.26 1.26 e, 0.9 0.9 TÂ„ (K) 2316 2316 2316 lf(J/g) 1046 25 1046 o (deg)* 135 135 h(W/cm2-K) 5 1.1x10-^ h (cm) 0.02 0.02 AT(K) 20 20 Ta (K) 2216 316 Up (cm/s) 0.13 0.062 (1/K) 3x10-^ 3x10-^ 7 (dyne/cm) 0.69 0.69 a (W/cm K") 5.67x10 '^ 5.67x10 '^ H (g/cm-s) 0.59 0.59 f(cm-/s) 0.19 0.19 ' (t>Â„ is defined as the angle between solid/melt interface and melt/gas interface.

PAGE 151

134 Table 4-3 Definition and magnitude of key dimensioniess parameters based on values given in Table 4-1 (Case 1: St = 1 and Case 2: St = 0.024). Case 1: St = 1 Case 2: St = 0.024 a.t Fo = Â— : Fourier number 2 0(1) 0(1) St = Â— ; Stefan number I J 1 0.024 Bo = ; Bond number 2.8x10^ 2.8x10-* Ra = ; Rayleigh number v^ 4.2x10"^ 4.2x10'^ 1 1 Arr, Ma = ; Marangoni number 1.56 1.56 Pe = ; Peclet number Â«/ 0.1 0.05 hr. Bi = Â— Â— ; Biot number k 1 with 0^ = -5 2.2x10-^ with 0^ = -100 Rod = ; Radiation number . K 0 1.8x10' ==^=-. ^

PAGE 152

135 Table 4-4 Two different scaling procedure and resulting nondimensional energy equations (a) choice 1 length scale: 1, = r^ , velocity scale: = a, / 1^ , time scale: t^ = Ir / , and temperature scale: AT = difference between temperature at melting point and at die tip. ae. dd. 1 a ae, a ae, energy eq, for melt Â— + F Â— = -Â—(RÂ—i) + Â— (Â— ) dx 'dZ RdR dR dZ dZ ae, ae, 13 ae, a ae, ensery for solid ^ * t/,^ [||(Â«^) ^ ^/ D a/f interface movement eg. Â— -k Â— = -t(1/ Â— )Â„ ^ dN dN St ^ dx ^ where R = r/1, , Z = z/1, , Up = Up/v, , V, = v,/v, , r = t/ t, , 0, = (T.-TJ/AT, = (T,-TJ/AT, p = p, / p, , and a = a, / a, (b) choice 2 length scale: 1, = rb , velocity scale: v, = St a, / 1, , time scale: t, = 1^ / v^ , and temperature scale: AT . at ^az /^a/? br bz dz energy eq. for solid St(^ + F-^) = a [1-^{R^) + Â— (-^)] ax 'dZ RdR dR dZ dZ interface movement eq. Â— k Â— = o (U Â— aiv aiv p dx ^

PAGE 153

136 Table 4-4 continued Boundary conditions based on choice 1 (1) for the melt: dQ, symmetric b.c. : Â— -=Q @R = Q, OiZ iH(Q) dR inlet temperature : 6, = 1 , @ 0 i i/?^ , Z = 0 80 jT T ' convectivelrcu&ative b.c. : " ^ = Bif,e, + ^^(6,+-^* (e^*-^*] , @ R = FiZ) , 0 i Z nH^ (2) for the solid: dd symmetric b.c. : Â— =0 , @ R = 0, //(O)
PAGE 154

0 137 meit 1 solid \^ j % =0 S(X,t) Governing equations: energy eq. for melt dd _ d^Q ^ 5t dX^ dY^ _ 50 ^ p^dS dN ~ St dx interface movement eq. Â— = -iÂ— B.C. 8 = 1, X=Q, T 6=0, XnS, t I.e. e(x,t =0.1) = 1 er/[X) 5(X.T,=0.1) = 2Ay^ Neumann's solution: 6 = 1erf^X) S = 2X^ ke^'erf{,X) = ^ Figure 4-1 Exact solution for melting of a solid confined region. to a semi-infinite

PAGE 155

138 Figure 4-2 Comparison of interface trajectories predicted by the coupled (full equations) and decoupled (quasi-stationary equations) approaches with three values of the Stefan numbers.

PAGE 156

139 Figure 4-3 Comparison of computed and exact solutions for St = 0.1303 St = 1.2216 and St = 2.8576. The figures on the left are the ' exact and computed temperature fields at different time instants. The figures on the right show the superposed exact and computed interface locations versus time. 1

PAGE 157

140 (a) schematic configuration in the melt region Qs Cp, U, T, k,^ = Gs (b) boundary condtions Figure 4-4 Schematic configurations and boundary conditions of tiie present EFG process.

PAGE 158

141 (a) Ap = 5 % 1^ 0.73 0.8 0.83 0.9 0.93 1 f/rb (b) Ap = 0 % Figure 4-5 Comparison of linear meniscus profile with the LaplaceYoung solution with R, = I, H, = 0.188. = 0.763, and various values of Bond number, (a) Ap = 5 % and (b) Ap = 0 %.

PAGE 159

142 (a) grid system Figure 4-6 J^;^^ ^^i^.,.r:n.^ coÂ„.ou. of U,e inicia, coÂ„di.ioÂ„s

PAGE 160

143 (a) grid system (b) isotherm (A0 = 0.2) Figure 4-7 The grid system and isothermal contours of the steadv-state solution for Case 1: St = 1. '

PAGE 161

144 (a) grid system (b) isotherm (A0 = 0.2) Figure 4-8 The grid system and isothermal contours of the steady-state solution for Case 2; St = 0.024.

PAGE 162

145 a e CO O II z > CJ C u O o 00 ._ h. o U f > o e O j= o c II o ^ CO Â£ Â•Â• a u .S II Â— Â«J O U O u O CO o W O I 3 OP

PAGE 163

CHAPTER 5 DYNAMIC SIMULATION OF THE EFG PROCESS 5.1 Introduction In the EFG technique, as in other meniscus-defined growth methods, the processes of heat and mass transfer serve as the link between the crystal growth condition and its resulting shape and structure. For example, the radius of the crystal is controlled by the interaction of heat transfer from the melt to the crystal and the ambient, with the shape of melt/gas meniscus controlled by the force balance between surface tension, the pressure difference across the interface, and the local shear stress distribution. Accordingly, the radius of the crystal can be affected by either the rate of heat transfer from the system to the surroundings, by varying the temperature of the melt that enters the die, or by changing the crystal pulling rate. In practice, both the die tip temperature and the puller speed may experience either intentional or unintentional variations in time. The unintentional variations can result from the perturbations to the operating conditions caused by environmental fluctuations, coil power fluctuations, or motor speed irregularities. The intentional variations can be designed to compensate the aforementioned perturbations to achieve a controlled growth process and uniform crystal properties and dimensions. In order to minimize the variations in fibre 146

PAGE 164

147 quality caused by these fluctuations in manufacturing conditions some type of active control scheme may be necessary (Gevelber et al. 1988, 1987, Gevelber and Stephanopoulos 1987, Meric 1979 and Robinson 1971). Furthermore, so far in the literature, the static contact condition is normally adopted. Effort has been made to develop two different models for the contact condition at the solid/melt/gas trijunction point; they are: (1) the conventional static model which requires that the Young's equilibrium contact condition, namely, a fixed contact angle between the surfaces separating melt/gas and solid/melt be maintained (Crowley 1983, Ettouney and Brown 1983, and Dussan V. 1979), and (2) a dynamic model which, instead of fixing the contact angle, allows it to vary according to the magnitude and direction of the instantaneous movement of the trijunction point. The purpose of this chapter is to mathematically analyze the transport phenomena in the EFG system so as to establish the connections between the actual growth conditions and the corresponding solidification characteristics. Two different values of Stefan number, St = 1 and 0.024, have been computed for the simulation. For St = 1, both temperature fluctuations at the die tip and speed fluctuations of the puller motor are considered; for St = 0.024, the pull speed fluctuations are investigated. Also addressed in the St = 0.024 case is the behavior of the trijunction point constrained by a dynamic contact condition. Relevant issues identified and insight gained from the simulation are discussed.

PAGE 165

148 5.2 Dynamic Perturbation on the EFG Process with St = 1 The steady-state solution for St = 1 shown in Figure 4-7 is used as the initial condition for the dynamic simulations. Two types of external perturbations will be considered: the fluctuation of temperature and, accordingly, the incoming heat flux, at the lower boundary of the melt, and the fluctuation of pull speed. The fluctuations will be represented either by combinations of sinusoidal waves or by a saw-toothed profile. 5.2.1 Base Temperature Perturbation We first consider the dynamic characteristics of the present EFG process in response to a base temperature perturbation. Varying the inlet temperature changes the amount of heat flux into the melt which, in turn, affects the solidification rate, the meniscus height H(R), and through the contact angle constraint, the trij unction location R,. In this section we examine the effect of the base temperature perturbation by keeping the pull speed (Up = 0.1,Pe =0.1) and all other operating parameters unchanged. The disturbance of the base temperature is given as ^(^) =(U.,-i;^..sin(^) (5.1) .-1 Wf where (0b)avg is the average base temperature , i.e. (6^),^^ = 1, Aj is the amplitude of the disturbance of the i* harmonic, and Ui is the period of the disturbance of the i"" harmonic.

PAGE 166

149 A time step size of At = 10'^ is used for all the calculations. Two forms of perturbation will be considered here; one with a single harmonic and another with three harmonics. For the perturbation with a single harmonic, the period, fi, is either 500 At or 20 At, both with the same amplitude, A =0.1. For the perturbation with three harmonics, the periods of the three harmonics are, respectively 500 At , 292 At and 108 At, each with the same magnitude of 0.1/3. In both cases, as shown in Figures 5-1 and 5-2, clear correlations exist between . He and R,. Firstly, while H, closely follows the movement of d^, there is a phase angle between them. The phase angle is 20 At (14.4Â°) for a single harmonic perturbation of the period of 500 At, 8 At (144Â°) for a single harmonic of the period of 20 At, and 12 At (8.6Â°) for the case of three-harmonic perturbation. The phase angle between and 6^, indicates the time required for the temperature perturbation at the base of the meniscus to travel to the melt/solid interface. Figure 5-1 also reveals that for both cases, and Rj are out of phase. This feature is a direct outcome of the constraint imposed at the trij unction point. As indicated by Eq. (4.41), with a constant Up and a fixed contact angle between the solid/gas and melt/gas interfaces, the rate of change of and R^ are of the opposite signs. Another aspect worth mentioning is the noniinearity of the system. Although both H, and R, move at the identical frequencies as those of their time-averaged mean values are different from the steady-state values. This phenomenon indicates that there are accumulated effects caused by the nonlinear

PAGE 167

150 physics intrinsic to the process. Finally, Figure 5-1 also clearly demonstrates that there is a large difference in the system's sensitivity, in terms of the magnitudes of He and Re variations, to low and high frequency forcings. The interface shapes and locations of the melt/solid at selected time instants are plotted in Figure 5-3 for both the single-harmonic and the threeharmonic forcings. While the maximum magnitudes of the temperature fluctuations in both cases are the same, the ranges of the solid/melt interface movement are not identical. The case with a single-harmonic perturbation (Q=500At) exhibits more interface movement than that with three-harmonic perturbations. This phenomenon largely reflects the fact that each harmonic affects the interface movement in a different time scale. It should be noted that due to the multi-dimensional nature of the problem, the heat flux difference across the solid/melt interface is not uniform, resulting in varying interface speeds, as evidenced by the nonuniform distances between the interfaces shown at different time instants shown in Figure 5-3. To better understand the effect of the forcing frequency of 6^, on the solid/melt interface movement, the phase angles among He , Up , and are inspected. By keeping the amplitude of the size of fluctuation unchanged, i.e. A =0.1, and varying the period from 500 At to 20 At, the computed results are summarized in Figure 5-4. Figure 5-4 (a) shows that the phase shift between He and Up increases as the frequency of the disturbance increases. He and Re are always out of phase because the rates of change of He and Re are of the opposite

PAGE 168

151 signs when Up is fixed, as the trijunction constraint, Eq, (4.41), implies. The small phase angle variation between and is caused by the solid/ melt interface curvature at the trijunction point. The percentages variation of H,. and with respect to the forcing frequency are presented in Figure 5-4 (b). Both the variations of H, and decrease as the forcing frequency increases. However, the characteristics of and R^ variations behave differently with respect to the forcing frequencies. For example, for the time period of &b perturbation larger than 200 At, the solidification process responds to the external forcing with approximately the same sensitivity. As the forcing frequency becomes sufficiently high, however, the sensitivities of both H, and R, decrease noticeably, indicating that the thermal inertia causes the system not to be able to keep pace with the high forcing frequencies. 5.2.2 Pull Speed Perturbation The pull speed in the typical EFG process (2-8 cm/min) is one to two orders of magnitude larger than that of conventional growth methods, e.g. Czochraiski (Zulehner 1983) or floating-zone (Duranceau and Brown 1986, and Ostrach 1980, 1983). For the process considered here, the radius of a grown crystal typically is only a few microns. As already discussed, convection is not important for the thermal transport in either the melt or the crystal. However, variation in pull speed has a direct impact on the solid/melt interface and the trijunction point movements, as reflected by Eq. (4.41). Increasing pull speed

PAGE 169

152 affects both the amount of heat transported by advection and the amount of latent heat released upon solidification, which cause the meniscus height to increase and the solid/melt interface to deform. For the results to be presented in the following, \]^{t) is perturbed while all other parameters are unchanged. The steady-state condition represented in Figure 4-7 is used as the initial condition to start the computations. Two types of perturbation are considered, namely, the sinusoidal and the saw-toothed shapes. With the sinusoidal perturbations, Up(T) varies as where (Up),vg = 0. 1 is average pull speed. This choice of the shape of perturbations is motivated by the desire to investigate the system dynamics with respect to (1) few discrete harmonics, and (2) a collection of harmonics with a wide distribution of frequencies. We first study two cases of a single-harmonic perturbation with J] = 500 At and 20 At, respectively, both with the same fluctuating magnitude, A = 0.5. The time histories of Up , H, and are shown in Figure 5-5. While both H, and R, respond to the Up perturbation at the same frequency, which is similar to the 6^, forcing, the decrease in the system's sensitivity with respect to forcing frequency (5.2) and with the saw-toothed perturbation, Up (t) varies as U,ir) = (UX, rim) forO
PAGE 170

153 seems more pronounced for U, than for 6, . There is also a clearly identifiable phase angle between H, and , which for the period of 500 Ar is 45 Ar (32.4^) and for 20 Ar is 5 Ar (90Â°). This phase angle represents the time scale needed by the thermal field in the melt to respond to the perturbation imposed by ; it is larger than that caused by the base temperature fluctuations with a low frequency perturbation, but smaller with a high frequency perturbation. The other noticeable aspect is that H, and R, are no longer out of phase, which is different from the case of d, perturbations. Again, the reason for this phenomenon can be found from the constraint at the trijunction point. Previously, Up was a constant, resulting in the requirement that dH,/dr and dR,/dr are essentially of opposite signs. Here, Up is no longer unchanged, and accordingly, the relationship between dH,/dr and dR,/dr is no longer unchanged. The phase angle between the time histories of H, and is a reflection of this characteristic. By keeping the forcing amplitude unchanged, i.e., A = 0.5, and varying the period from 500 Ar to 20 Ar, the relative phase angles between H, and Up, between H, and K , as well as the percentage variations of H, and RÂ„ are depicted in Figure 5-6. Figure 5-6(a) shows that the phase shift between H, and Up increases as the forcing frequency increases. The phase shifts between both H, and Up , and between H, and R, become larger as the frequency becomes higher. Figure 5-6(a) indicates that H, and R, tend to be more out of phase as the forcing frequency increases. The trend of the phase angle between H, and external forcing is largely the same, as evidenced by the similarity between -Figures 5-4(a) and 5-6(a). The variations of

PAGE 171

154 H, and with respect to the forcing frequency are presented in Figure 5-6(b). Sensitivities of both H, and R, decrease as the forcing frequency increases. It indicates that by perturbing either 6^, or Up, the thermal inertia of the system cannot respond equally with respect to different harmonics. However, unlike the case of By, forcing shown in Figure 5-3(b), Figure 5-6(b) indicates that the sensitivities of both H, and R, reduce more or less at the same rates as the forcing frequency of Up increases. We next study the movement of the solid/melt interface in response to the perturbed Up containing three incommensurate harmonics, with the periods of 500 At, 292 At and 108 At, all with the same magnitude of 0.5/3. The frequencies imposed here are identical to a case of perturbing 6^, studied previously, shown in Figure 5-2. Figure 5-7 shows the time histories of Up , H, and R, for this case. While Up in Figure 5-7 and 6^ in Figure 5-2 are of the identical form, the shapes and magnitudes of variations in H, and R, are quite different in the two cases. With the 0b perturbation, all three forcing frequencies are better captured by both He and R^ . As already observed, by perturbing Up , both H, and R, progressively reduce their sensitivities with respect to the external forcing; by penurbing ^^Â„ on the other hand, H,, and especially RÂ„ maintain their sensitivities with respect to the forcing frequency until it becomes quite high. To further illustrate this characteristic, a saw-toothed perturbation is given to both Up and 6^, with a period of 20 At and A =0.1. The time histories of Up, He and R^ are shown in Figure 5-8. With such high frequency perturbations, while

PAGE 172

155 the system responds to the forcing in a modest manner, the sensitivity to Up is much smaller than to Â• 5.3 Dynamic Perturbation on the EFG Process with St = 0.024 Next, we choose a Stefan number of 0.024, which is more representative of the current state-of-the-art in actual production. Together with previous section St = l,the dynamic characteristics of the EFG process under different conditions can be assessed and compared. For the results presented here, all the relevant conditions are given in Table 4-3. Compared to the case of St = 1, the steadystate pull speed is reduced by half and some modification are made to the thermal boundary conditions. The motivation for making these changes is to bring the modeled condition closer to a typical production condition. The most important change is the Stefan number, other modification do not change the qualitative features presented below. The simulation starts with the steady-state solution. Figure 5-8, as the initial condition. Figure 5-9 shows the time series plots of two cases, one with a fluctuating pull speed of the period of 500At, At=10'' ,and another with 20 At. Shown there are the time histories of the pull speed. Up , the height and the radius of the trijunction point, H, and R, , respectively. Consistent with the case of St = 1 studied in previous section, the crystal diameter responds to the externally imposed perturbation at the corresponding frequency, but the sensitivity of the response depends on the frequency of perturbation. For a slower perturbation of the period of 500At, both

PAGE 173

156 H, and R, exhibit higher levels of oscillation than for a faster perturbation of 20At; the differences between them are substantial. Furthermore, also consistent with the case of St = I, the nonlinearity of the underlining physics has caused the time averaged values of H, and R, to be different from the steady-stale values. Figure 5-10 gives an overall depiction of the percentage variations of H, and R, with respect to different perturbation periods predicted by the present model. For both and R, , the tendency of exhibiting reduced sensitivities to the external perturbations as the frequency increases can be clearly observed. The reason for this phenomenon obviously lies in the relative competition between the time scale of the grower and the time scale of the perturbation. Since At is 10'^ , for the pull speed fluctuation of the period of several hundred At, the system's dominant non-dimensional time scale, which is of order 1 according to choice 2 given in Table 4-3, and the perturbations time scale are comparable. Hence, the system is able to respond and follow the perturbation closely, resulting in a clear connection between the dynamic behaviors of Up and H, and R^ . For high frequency perturbations, of the period of hundred At or less, the system's time scale is too long, causing decreases in sensitivity of and Rj , as observed in Figure 5-10. This phenomenon has also been observed in the actual growth system. Figure 5-11 shows a Bode diagram based on a series of experiments to depict the frequency response between the fibre diameter and pull speed (Backman 1992). It indicates that with relatively slow varying perturbations over a range of frequencies in pull speed, the fibre diameter responds with a

PAGE 174

157 certain sensitivity; as the perturbation frequency increases beyond a threshold, the system can no longer respond to it with equal sensitivities. Hence, the predicted trend of the system response is qualitatively verified experimentally. 5.4 Dynamic Contact Angle While it seems reasonable to use the LaplaceYoung equation to determine the meniscus shape due to the smallness of Ra, Pe and Ma, it is an appropriate prescription of the contact condition at trijunction point. Two models have been investigated in this regard, including a conventionally adopted static model which fixes the contact angle, and a dynamic model which regulates the contact angle according to the speed and direction of the instantaneous movement of the trijunction point. The primary focus of this section is to investigate the effect of the perturbation in pull speed, subject to two different models for the contact condition at trijunction point, on the characteristics of the dimensional variation of the fibre. We present the results obtained based on a dynamic model for the contact condition at the trijunction point. The model, as schematically depicted in Figure 5-12, now allows different values of contact angle according to the instantaneous direction as well as the speed of movement of the trijunction point. Drawing the analogy to the advancing and receding contact angle of a liquid drop on a flat surface (Dussan V. 1979), two different contact angles are assigned as the asymptotic values for the outward and inward moving cases, respectively.

PAGE 175

158 Furthermore, if the contact point moves at a slow enough speed, in the range of +e and -e as shown in Figure 5-12, then the contact angle is allowed to vary between the two asymptotic angles. For the case implemented, we take the asymptotic value of the contact angle associated with the outward moving case to be the same as that adopted in the static model, namely, 135Â°. The asymptotic contact angle of the inward moving case is taken as 134Â°, and e is taken as 0.1. In other words, the variation of contact angle between the inward and outward moving cases is 1Â°. While this variation seems very small, it has a quite noticeable impact on the characteristics of the solidification process. Figure 5-13 shows the comparison of the time histories of both H, and between the cases with static and dynamic contact conditions, where a perturbation of Up with a period of 100 At is enforced. It is striking to observe that with a seemingly small modification of the contact condition, the resulting values of H, and R, are substantially affected. For R, , both time-averaged and fluctuating values are quite different with the two contact conditions; the dynamic contact model yields a crystal of a somewhat larger mean diameter and smaller variations. For , on the other hand, the time-averaged fibre diameter is smaller with the use of the dynamic contact condition; the fluctuating magnitudes in both cases, however, appear quite comparable. Hence, there seems to be differences in behavior between H, and R, resulting from the static and dynamic models. In this regard, several points can be noted to elucidate the results observed here. As observed in Figure 4-8, the solid/melt interface is concave

PAGE 176

159 toward the solid phase, indicating that as the contact angle is reduced, the trijunction point tends to move outward, yielding a larger diameter as shown in Figure 5-13. This will also seem to cause a corresponding increase in H, , as seen in the initial transient behavior in the time series plot of . However, a larger crystal diameter also causes a larger volume and a higher heat loss both via conduction and advection. Accordingly, there is also a change in the balance of the thermal constraint at the interface, which eventually reduces the value of H, as observed in the later stage of Figure 5-13. It should be noted that if the solid/melt interface were of the opposite curvature, then the trend of R, will change as well. Hence, the interaction of thermal and contact condition plays a pivotal role in determining the solidification characteristics. Furthermore, since the contact angle is no longer fixed in the dynamic model, an extra degree of nonlinearity also appears via the enforcement of the contact condition. Plots of the rate of change of H, and R, are given in Figure 5-14 to demonstrate this point. Clearly, with the static contact condition, a single harmonic reflecting that of the imposed perturbation of Up is reproduced in both H, and R, . On the other hand, with the dynamic contact condition, the movement of trijunction point is no longer consistent with the imposed frequency, instead, it exhibits a quasi-periodic pattern with many additional frequencies present. These additional frequencies appear because the way the contact angle, and hence the values of H, and R, , vary in response to the pull speed perturbations.

PAGE 177

160 In general, the value of the equilibrium contact angle may be difficult to quantify. Here an assessment is made to investigate the effect of contact angle on the solidification dynamics. Three values of static contact angle, namely, 107Â°, 125Â° and 135Â°, are used. Figure 5-15 shows the frequency responses of solid/melt interface with different static contact angles, to the pull speed fluctuations. The overall features of the interface movement of the three contact angles show similar trend for both AH, and AR, , although the variations of AH, and AR, are affected by the value of the contact angle. Because the curvature of the interface near the trij unction point is sensitive to the local energy balance and does not vary monotonically with the contact angle, there is no simple trend between AH, , ARj and
PAGE 178

161 frequencies. It has also been predicted that as the forcing frequency increases, the solid/melt interface responds at a reduced sensitivity. This trend has been observed experimentally. The time-averaged values of and are different from the steady-state values because of the nonlinearity of the fibre growth process. 2. While the Stefan number substantially affects the temperature distribution, interface locations, and crystal diameter, it does not affect the qualitative behavior of the solidification process in the dynamic environment. 3. The phase shifts as well as magnitudes of H,. and are found to be forcing frequency dependent for both base temperature and pull speed perturbations, (i) The phase angles between H, and external perturbation increase as the forcing frequency increases, (ii) Both H, and R, are less affected by the high frequency harmonics of the perturbations. However, with the Up perturbation, responses in H,. and R^ decrease at relatively uniform rates as the forcing frequency increases, but with the perturbation, and R, exhibit two different sensitivity scales with respect to the forcing frequency. 4. The relationship between H, and R, is dictated by the constraint at the trijunction point, which differs in two different types of forcing mechanisms studied. With the base temperature perturbation, while and R, remains essentially out of phase at all frequencies. H, and maintain a phase angle which is frequency dependent. With the pull speed penurbation, the relationship between He and R, , as well as between H^ and Up , is frequency dependent.

PAGE 179

162 5. Two different models have been utilized for the contact condition at trij unction point, including a conventionally adopted static contact angle, and a dynamic model which allows the contact angle to vary according to the direction and the rate of the trij unction point movement. It turns out that the two models can yield substantial differences in the crystal diameter, including the timeaveraged as well as the fluctuating magnitudes. Furthermore, additional harmonics can also be yielded by the dynamic contact condition due to the extra nonlinearity contained by it. These observations have clearly indicated that the static contact condition is not sufficient to yield all the information regarding the dynamic behavior of the crystal correctly. 6. Available experimental information has confirmed the characteristics of frequency response of discussed here.

PAGE 180

163 (a) Q = 500 At 0 210 400 Â«00 500 1000 1200 I40O 1000 1100 2000 time steps (b) Q = 20 At to 100 120 140 leo lis 200 time steps time steps ao 10 ion 120 |40 loo no inn time steps 0.79 0.7Â» 0.77 0.7Â« ' 0.73 0.74 0.71 0.72, 20 Â« Â« io 100 iio 155 iÂ«o Â— lao ino time steps Figure 5-1 Time histories of H, and R, subject to a single harmonic perturbation of with (a) Q = 500 Ar and (b) Q = 20 At (St = 1.0)

PAGE 181

Figure 5-2 Time histories of and R, subject to a three-harmonic perturbation of with 12, = 500 Ar, ^ = 292 At, and = 108 Ar . (St = 1.0)

PAGE 182

165 0^5 p 0.24 0^ 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 r (a) a singie-hannonic permrbaiion in 6^ with Qj = 500 At 0.21 tÂ»l300dt \. t1600dt Â• tÂ» 1700 dt "Â•Â• Â• ... i-1800dt n/'S... 0 0.1 Oa OJ 0.4 OJ 0.6 0.7 0.8 (b) a three-harmonic pemirbauon in 6^ with Q, = 500 At, ~ 292 At and Q, = 108 At . Figure 5-3 The interface shapes at different time instants of (a) a singleharmonic perturbation in d^, with Q, = 500 Atand (b) a threeharmonic perturbation in 6^, with Q, = 500 Ar, = 292 Ar, and = 108 Ar . (St = 1.0)

PAGE 183

166 130 -200 i , , . . I Â« 50 100 130 200 230 300 330 400 430 300 penod (lime stept) (Â«) Hie relative pluue uiBiei between He Â«Bdefc and between and R, 18 Ol ' < ^ , , , 0 30 100 130 200 230 300 330 400 430 300 period (time steps) (b) Petcentage variaiioni of andR^ Characteristics of interface movement with respect to different forcmg frequencies of d, (a) The relative phase angles between H and 5, and between, H, and . (b) Percentage variations of H, and R, . (St = 1.0)

PAGE 184

Figure 5-5 Time histories H. and subject to a single harmonic (St = ^^'^ " ^ """^ 0 = 20 Ar

PAGE 185

168 130 100 30 -30 -100 -130 -200 slaiive phase angle between AU^ pelatiye phase angle between 10 100 130 200 230 300 330 400 430 300 period (time steps) (at -n^ reladve phase angle, between H ^ and Up and between and R , 3J 2J IJ 0.3 % variation of 30 100 130 200 230 300 330 400 430 300 period (time steps) (b)Petcentago variations of He andRg Figure 5-6 Characteristics of interface -movement with respect to different H ani uTT: ^= "'^''^ P*-' between o/ tan
PAGE 186

169 ure 5-7 Time histories H, and subject to a three-harmonic perturbation of with fi, = 500 A., = 292 Ar and "3 108 At . (St = 1.0)

PAGE 187

170 (a) base temperanire pemirbatiaa (b) puU speed pematoaucm to 100 120 140 |Â«0 m n tt *" * " iw ijo MO ieo itS Â— no 10 vn o.n o.n 0.T7 0.7S . O.TSi 0.74 0L7J 0.7J 20 40 SO 10 100 120 140 IM \to' Figure 5-8 J^^^ories of saw-toot^^^ forcing with 1] = 20 A. and \ corresponding H, and R,. (St = 1 0) (a base temperature perturbation 6,(7), ^ ^ (b) pull speed perturbation U (r).

PAGE 188

171 (a) period = 500 Ax (b) period = 20 At rime steps time steps time steps Figure 5-9 Time histories of and subject to a single harmonic perturbation of U, with (a) period=500 Ar and (b) period =20 Ar. (St = 0.024)

PAGE 189

172 c u 1 ; I : I I : 3 U op ^ pas JO QOifBUBA aSesoaojad ^

PAGE 190

173 u 'o c a 3 U 60 Â— .E > c Â•Â— _ ^ CO a. Â« a ^ E D < -2 2 3 ^ m aopcuBA

PAGE 191

174 c o + to a CO u CO C o c > es ISO C n o 2 c o U CO 3 O c c o3 C S es O Â« c o u s '-B u o o. c o c -s.-* Â•s.'' -s? c o o u u e o o e 00

PAGE 192

175 0355 time steps 0.44 time steps Figure 5-13 Time histories of H, and subject to a single harmonic perturbation of Up with period = lOOAr . (St = 0.024)

PAGE 193

176 J I 1 ' 1 1 > 1800 1S20 1840 1860 1880 1900 1920 1940 1960 1980 ZOOO time steps 1800 1820 1S40 1860 1880 1900 1920 1940 1960 1980 2000 time steps ure 5-14 Time histories of dH./dr and dR./dr subject to a single harmonic perturbation of Up with period = lOOAr . (St = 0.024)

PAGE 194

177 o c o > u 00 B c u u CU c Â•c > u OA iS c u u 0 50 100 150 200 250 300 350 400 450 500 Period (non-dimensional time steps) (a) Percentage variation of 12 10 8 4 50 100 150 200 250 300 350 400 450 500 Period (non-dimensional time steps) (b) Percentage variation of R, Figure 5-15 Characteristics of interface movement with respect to different forcing periods of Up for different static contact angles (St = 0.024)

PAGE 195

CHAPTER 6 CONCLUSIONS AND RECOMMENDATIONS This chapter presents some conclusions that have been drawn from this work, together with a few suggestions as to how this work could be usefully extended. 6.1 Achievements and Findings An adequate theoretical formulation of the EFG process needs to address several issues. They include the method of predicting the thermal fields in the solid and the liquid, the shapes of solid/liquid and liquid/gas interfaces, front tracking and moving grid techniques, understanding of the meniscus behaviors and the contact condition at the trij unction point, and appropriate scaling procedures. In this study, a thermocapillary model has been developed to investigate the dynamic characteristics and associated transport phenomena of the EFG process, subject to the external perturbations and different contact conditions. Based on the computed results, the following conclusions are made. An order of magnitude estimation shows that the convective mechanisms are minor influences on the transport process in a typical EFG process; conduction dominates energy transport. In the absence of an externally imposed pressure field, capillarity dominates the shape of the meniscus. 178

PAGE 196

179 The shape of the liquid/gas meniscus that connects the growing crystal with the die tip, governed by the LaplaceYoung equation and subject to the contact condition at the trijunction point, is critical in determining the properties of the solidified crystal. A study of meniscus formation shows that multiple solutions may exist. The realizability of these possible meniscus profiles has been examined using the minimum free energy concept. Interactions among Bond number, pressurization, aspect ratio, contact condition, and pulling direction have been found to result in a variety of meniscus behaviors. A hysteresis model for the dynamic contact line motion has been employed to investigate the quasiequilibrium dynamics of the meniscus. It is found that the relation in terms of waveform and phase of the radius of the trijunction point with respect to the height of the trijunction point cannot be fully captured by the formula conventionally employed in crystal growth simulations. Full fluid-dynamic consideration is required to assess the implication of such behavior. A new generalized interface tracking technique is developed, which overcomes restrictions of single-valuedness of the interface imposed by commonly used mapping methods. This method is shown to be both accurate and versatile. Combining this method with a moving grid technique, detailed thermal fields in both solid and liquid phases can be computed. Two different scaling procedures for the interface movement are employed for both St = 1 and St = 0.024. One scaling procedure uses the diffusional characteristic velocity as the velocity scale, and the other includes both the

PAGE 197

180 diffusional velocity and the latent heat effects. Since the interface movement is primarily responsible for the nonlinearity of the phase change problem, the procedure which accounts for the latent heat effect results in better computational efficiency for low St cases. The sensitivities and dynamic characteristics of the EFG system subject to the base temperature and pull speed perturbations are explored. For both disturbances, the crystal radius and the height of the trijunction point H, reach limit cycles at the same frequency as the external forcing; the magnitude of their variations with respect to these perturbations are found to be dependent on the forcing frequencies. As the forcing frequency increases, the solid/liquid interface responds with a reduced sensitivity. This trend has been observed experimentally. Also, the time-averaged values of and R, are different from the steady-state values because of the nonlinearity of the solidification process. Two different models have been utilized to simulate the contact condition at the trijunction point, including a conventionally adopted static contact condition, and a dynamic model in which the contact angle varies according to the direction and the simultaneous speed of the trijunction movement. It turns out that the two contact conditions can yield substantial differences in solidification characteristics as reflected by the variation in crystal diameter. These results indicate that the conventional static model is not adequate to model some aspects of the trijunction point behavior, particularly in the presence of the phase change.

PAGE 198

181 6.2 Directions for Future Research The present thermocapillary model has produced some interesting results, and has been partially confirmed experimentally. It will offer useful insights into a practically important process. To further extend the content of the model, further work can be conducted in research areas as suggested below. The contact condition specified at the trij unction point is found to be important in determining the crystal growth characteristics. Our simplified dynamic contact model is a first step toward accounting for the nonequilibrium effects at the trij unction point. Hydrodynamic problems involving dynamic contact lines are unsteady, highly nonlinear, and involve free and moving interfaces with singular solution near the moving contact lines (Shopov and Bazhlekov 1991, Shopov et al. 1990, Dussan V. 1979). An adequate theory pertaining to such systems is still in the formative stages, and much work remains to be done. Experimental information is critical in helping to gain a better understanding of this issue. In this study, the phase boundary is determined by the condition of local thermodynamic equilibrium. This assumption is only valid for small undercooling and slow pull speed. In general, the interface temperature is a function of the local rate of phase change kinetics, introducing an extra degree of freedom which is constrained by a kinetic relationship between interface velocity and undercooling (Wang and Matthys 1992, Clyne 1984, Levi and Mehrabian 1982, Shingu and Ozaki, 1975). A complete model capable of predicting accurately the

PAGE 199

182 solid/liquid interface velocity and interface location according to both the thermal and kinetic constraints will be useful. On a practical side, the findings reported here can be used to facilitate design of appropriate active control schemes to improve the fibre quality (Gevelber et al. 1988, 1987, Gevelber and Stephanopoulos 1987, Meric 1979, Robison 1971).

PAGE 200

REFERENCES Albert, M. R. and O'Neill, K., "Moving Boundary-Moving Mesh Analysis of Phase Change Using Finite Element Method," International Journal for Numerical Methods in Engineering, 23, pp. 591-607, 1986. Anderson, D. A., Tannehill, T. C. and Fletcher, R. H., Computational Fluid Mechanics and Heat Transfer . Hemisphere, Washington, D.C.,1984. Backman, D. G., "Intelligent Processing of MMC Materials," Final Report to DARPA/NRL, Contract No. N00014-90-C-0060, GE Aircraft Engines, Lynn, MA, 1992. Bankoff, S. G., "Heat Conduction or Diffusion with Change of Phase," in: Advances in Chemical Engineering . Drew, T. B., Hoopes, J. W. Jr., Vermeulen, T. and Cokelet, C. R. (eds.), 5, pp. 75-150, Academic Press, New York, 1964. Beckermann, C. and Viskanta, R., "Double-Diffusive Convection Due to Melting," International Journal of Heat Transfer and Mass Transfer, 31, pp. 2077-2089, 1988. Bennon, W.D. and Incropera, F. P., "A Continuum Method for Momentum, Heat and Species Transport in Binary Solid-Liquid Phase Change Systems-I. Model Formulation," International Journal of Heat and Mass Transfer, 30, pp. 2161-2170, 1987. Benton, E. R. and Platzman, G. W., "A Table of the Solutions of the OneDimensional Burgers Equation," Quarterly of Applied Mathematics, 30, pp. 195212, 1972. Boley, B. A. "The Embedding Technique in Melting and Solidification Problems," in: Moving Boundary Proble ms in Heat Flow and Diffusion . Ockendon J. R. and Hodgkins W. R. (eds.), pp. 150-172, Oxford University Press, London, 1975. Brent, A.D., Voller, V. R. and Reid, K. J., "Enthalpy-Porosity Technique for Modeling Convection-Diffusion Phase Change: Application to the Melting of a Pure Metal," Numerical Heat Transfer, 13, pp. 297-318, 1988. 183

PAGE 201

184 Brown, R. A., "Theory of Transport Processes in Single Crystal Growth from the Melt," AIChE Journal, 34, pp. 881-911, 1988. Carslaw, H.S. and Jaeger, J. C, Conduction of Heat in Solids . London, Oxford University Press, London, 1959. Chemousko, F. L., "Solution of Non-Linear Heat Conduction Problems in Media with Phase Change," International Chemical Engineering, 10, pp. 42-48, 1970. Clyne, T. W., "Numerical Treatment of Rapid Solidification," Metallurgical Transaction B 15B, pp. 369-381, 1984. Correa, S. M. and W. Shyy, "Computational Models and Methods for Continuous Gaseous Turbulent Combustion," Progress in Energy and Combustion Science, 13, pp. 249-292. 1987. Crank, J. . Free and Moving Boundary Problems . Oxford University Press, London, 1984. Crank, J. and Gupta, R. S., "Isotherm Migration Method in Two-Dimensions," International Journal of Heat and Mass Transfer, 18, pp. 1101-1106, 1975. Crochet, M. J.,Geyling, F. T. and Schaftingen, J. J. van, "Numerical Simulation of the Horizontal Bridgman Growth. Part 1: Two-Dimensional Flow," International Journal for Numerical Methods in Fluids, 7, pp. 27-49, 1987. Crowley, A. B., "Mathematical Modeling of Heat Flow in Czochralski Crystal Growth," IMA Journal of Applied Mathematics, 30, pp. 173-189, 1983. Cuvelier, C. and Schulkes, R. M. S. M., "Some Numerical Methods for the Computation of Capillary Free Boundaries Governed by the Navier-Stokes Equations," SIAM Review, 32, pp. 355-423, 1990. Dantzig, J. A., "Modeling Liquid-Solid Phase Change with Melt Convection," International Journal for Numerical Methods in Engineering, 28, pp. 1769-1785, 1989. de Gennes, P. G., "Wetting: Statics and Dynamics," Reviews of Modem Physics, 57, pp. 827-863. 1985. Derby, J. J. and Brown, R. A., "Thermal-Capillary Analysis of Czochralski and Liquid Encapsulated Czochralski Crystal Growth," Journal of Crystal Growth, 74, pp. 605-624, 1986.

PAGE 202

185 Dupont, S.,Marchal, J. M., Crochet, M. J. and Geyling, F. T., "Numerical Simulation of the Horizontal Bridgman Growth. Part 2: Three-Dimensional Flow," International Journal for Numerical Methods in Fluids, 7, pp. 49-67, 1987. Duranceau, J. D. and Brown, R.A., "Thermal-Capillary Analysis of Small-Scale Floating Zones: Steady-state Calculations," Journal of Crystal Growth, 75, pp. 367389, 1986. Dussan V.,E. B.. "On the Spreading of Liquid on Solid Surface: Static and Dynamic Contact Lines," Annual Review of Fluid Mechanics, 11, pp. 371-400, 1979. Dyson, D. C, "The Energy Principle in the Stability of Interfaces," Progress in Surface and Membrane Science, 12, pp. 479-564, 1978. Eick. J. D., Good, R. J. and Neumann, A. V., "Thermodynamics of Contact Angles. II. Rough Solid Surfaces," Journal of Colloid and Interface Science, 53, pp. 235248, 1975. Epstein, M. and Cheung, F. B., "Complex Freezing-Melting Interfaces in Fluid Flow, "Annual Review of Fluid Mechanics, 15, pp. 293-319, 1983. Ettouney H. M. and Brown, R. A., "Finite-Element Methods for Steady Solidification Problems," Journal of Computational Physics, 49, pp. 118-150,1983. Fletcher, C. A. J.. Computational Techniques for Fluid Dvnamics . Vols. I and II, Springer, New York, 1988. Floryan, J. M. and Rasmussen, H., "Numerical Methods for Viscous Flows with Moving Boundaries," Applied Mechanics Reviews, 42, pp. 323-340, 1989. Gau, C. and Viskanta, R., "Melting and Solidification of Pure Metal on a Vertical WaU," Journal of Heat Transfer, 108, pp. 174-181,1986. Gevelber, M. A. and Stephanopoulos, G., "Dynamics and Control of the Czochralski Process I. Modelling and Dynamic Characterization," Journal of Crystal Growth, 84, pp. 647-668, 1987. Gevelber, M. A., Stephanopoulos, G. and Wargo, M. J., "Dynamics and Control of the Czochralski Process II. Objectives and Control Structure Design," Journal of Crystal Growth, 91, pp. 199-217, 1988. Gevelber, M. A.. Wargo, M. A. and Stephanopoulos, G., "Advanced Control Design Considerations for the Czochralski Process, " Journal of Crystal Growth, 85, pp 256263, 1987.

PAGE 203

186 Goodman, T. R., "The Heat Balance Integral and Its Application to Problems Involving a Change of Phase," ASME Journal of Heat Transfer, 80, pp. 335-342, 1958. Haley, P. J. and Miskis, J. M., "The Effect of the Contact Line on Droplet Spreading," Journal of Fluid Mechanics, 223, pp. 57-81, 1991. Heinrich, J. C, Felicelli, S.,Nadapukar, P. and Poirier, D. R., "Thermosolutal Convection during Dendritic Solidification of Alloys, Part II: Nonlinear Convection," Metallurgical Transaction B 20B, pp. 883-891, 1989. Hirsch, C. Numerical Computation of Internal and External Flows . Vol. 2, Wiley, New York, 1990. Huh, C. and Mason, S. G., "Effects of Surface Roughness on Wetting," Journal of Colloid and Interface Science, 60, pp. 11-38, 1977. Huh, C.,and Scriven, L. E., "Shapes of Axisymmetric Fluid Interfaces of Unbounded Extent," Journal of Colloid and Interface Science, 30, pp. 323-331, 1960. Imaishi, N.,Tsukada, T.,Hozawa, M.,Okano, Y. and Hirata, A., "Numerical Study of the Czochralski Growth of Oxide Single Crystal," in: Heat and Mass Transfer in Material Process . Lior, N. and Tanasawa I. (eds.), pp. 123-136, Hemisphere, New York, 1991. Jansons, K. M., "Moving Contact Lines at Non-zero Capillary Numbers," Journal of Fluid Mechanics, 167, pp. 393-407, 1986. Joanny, J. F. and de Gennes, P. G., " A Model for Contact Angle Hysteresis," Journal of Chemical Physics, 81, pp. 552-562, 1984. Johnson, R. E., "Conflicts between Gibbsian Thermodynamics and Recent Treatments of Interfacial Energies in Solid-Liquid-Vapor Systems," Journal of Physical Chemistry, 63, pp. 1655-1659, 1959. Kalejs, J. P., Chin, L.-Y. and Carlson, F. M., "Interface Shape Studies for Silicon Ribbon Growth by the EFG Technique," Journal of Crystal Growth, 61, pp. 473484, 1983a. Kalejs, J. P.,Ettouney, H. M. and Brown, R. A., "Comparison of Growth Characteristics of Sapphire and Silicon Ribbon Produced by EFG," Journal of Crystal Growth, 65, pp. 316-32, 1983b.

PAGE 204

187 Katoh, K.,Fujita, H. and Sasaki, H., "Macroscopic Wetting Behavior and a Method for Measuring Contact Angles, " Journal of Fluids Engineering, 112, pp. 289-295, 1990. Koplik, J.,Banavar, J. and Willemsen, J. P., " Molecular Dynamics of Fluid Flow at Solid Surfaces," Physics of Fluids, A 1, pp. 781-794, 1989. LaBeile, H. E. Jr., "EFG, the Invention and Application to Sapphire Growth," Journal of Crystal Growth, 50, pp. 8-17, 1980. Lacroix, M., "Computation of Heat Transfer During Melting of a Pure Substance from an Isothermal Wall," Numerical Heat Transfer, B, 15, pp. 191-210, 1989. Lacroix, M. and VoUer, V. R., "Finite Difference Solutions of Solidification Phase Change Problems: Transformed versus Fixed Grids," Numerical Heat Transfer, B, 17, pp. 25-41, 1990. Leonard, B. P., "A Stable and Accurate Convective Modeling Procedure Based on Quadratic Upstream Interpolation," Computer Methods in Applied Mechanics and Engineering, 19, pp. 59-98, 1979. Levi, C. G. and Mehrabian, R., "Heat Flow during Rapid Solidification of Undercooled Metal Droplets," Metallurgical Transaction, A 13 A, pp. 221-234, 1982. Lightfoot, N. M. H., "The Solidification of Molten Steel," London Mathematical Society Proceedings, 31, pp. 97-116, 1929. Lobar, B. L. and Jain, P. C, "Variable Mesh Cubic Spline Technique for N-Wave Solution of Burgers' Equation," Journal of Computational Physics, 39, pp. 433442, 1981. Lynch, D. R., "Unified Approach to Simulation on Deforming Elements with Application to Phase Change Problem," Journal of Computational Physics, 47, pp. 387-441, 1982. Meric, R. A., "Finite Element and Conjugate Gradient Methods for a Nonlinear Optimal Heat Transfer Control Problem," International Journal for Numerical Methods in Engineering, 14, pp. 1851-1863, 1979. Michel, D. H., "Meniscus Stability," Annual Review of Fluid Mechanics, 13, pp. 189-215, 1981.

PAGE 205

188 Miller, C. A. and Neogi, P.. Interfacial Phenomena . Marcel Dekker, New York, 1985. Morton, K., "Stability and Convergence in Fluid Flow Problems," Proceedings of the Royal Society of London, Series A (Mathematical and Physical Science), 323, pp. 237-253, 1971. Neumann, A. V. and Good, R. J., "Thermodynamics of Contact Angles, I. Heterogeneous Solid Surfaces," Journal of Colloid and Interface Science, 38, pp. 341-358, 1972. Ngan, C. G., and Dussan V.,E. B., " On the Nature of the Dynamic Contact Angle: An Experimental Study," Journal of Fluid Mechanics, 118, pp. 27-40, 1982. Ockendon, J. R. and Hodgkins, W. R.. Moving Boundary Problems in Heat Flow and Diffusion . Oxford University Press, London, 1975. Oliver, J. F.,Huh, C. and Mason, S. G., "Resistance to Spreading of Liquids by Sharp Edges," Journal of Colloid and Interface Science, 59, pp. 568-581, 1977. Ostrach, S., "Fluid Mechanics in Crystal Growth 1982 Freeman Scholar Lecture," ASME Journal of Fluid Engineering, 105, pp. 5-20, 1983. Ostrach, S., "Low-Gravity Fluid Flows, " Annual Review of Fluid Mechanics, 14, pp. 313-345, 1980. Ozisik, M. N.. Heat Conduction . Wiley, New York, 1980. Ozisik, M. N., "A Note on the General Formulation of Phase Change Problems as Heat Conduction Problems with a Moving Heat Source, " ASME Journal of Heat Transfer, 100, pp. 370-377, 1978. Patankar, S. V.. Numerical Heat Transfer and Fluid Flow . Hemisphere. New York, 1980. Pelce, P.. Dynamics of Curved Fronts. Perspectives in Physics . Academic Press Inc., New York, 1988. Pitts, E., "The Stability of a Meniscus Joining a Vertical Rod to a Bath of Liquid," Journal of Fluid Mechanics, 76, pp. 641-651, 1976. Roache, P. J., "On Artificial Viscosity;" Journal of Computational Physics, 10, pp. 169-184, 1972.

PAGE 206

189 Robison, A. C, "A Survey of Optimal Control of Distributed-Parameter Systems," Automatics, 7, pp. 371-388, 1971. Roe, P. L., "Characteristic-Based Schemes for the Euler Equations," Annual Review of Fluid Mechanics, 18, pp. 337-365, 1986. Sachs, E. M., "Thermal Sensitivity and Stability of EFG Silicon Ribbon Growth," Journal of Crystal Growth, 50, pp. 102-113, 1980. Sackinger, P. A., Brown, R. A. and Derby, J. J., "A Finite Element Method for Analysis of Fluid Flow, Heat Transfer and Free Surface in Czochralski Crystal Growth," International Journal for Numerical Methods in Fluids, 9, pp. 453-492, 1989. Saitoh, T., "Numerical Method for Multi-Dimensional Freezing Problems in Arbitrary Domains," ASME Journal of Heat Transfer, 100, pp. 294-299, 1978. Salcudean, M. and Abdullah, Z., "On the Numerical Modelling of Heat Transfer during Solidification Process," International Journal for Numerical Methods in Engineering, 25, 1988, pp. 445-473. Schwartz, A. M. and Tejada, S. B., "Studies of Dynamic Contact Angles on Solids," Journal of Colloid and Interface Science, 38, pp. 359-375, 1982. Shingu, P. H. and Ozaki, "Solidification Rate in Rapid Conduction Cooling," Metallurgical Transaction A 6A, pp. 33-37, 1975. Shopov, P. J. and Bazhlekov, I. B., "Numerical Method for Viscous Hydrodynamic Problems with Dynamic Contact Lines, " Computer Methods in Applied Mechanics and Engineering, 91, pp. 1157-1174,1991. Shopov, P. J., Minev, P. D., Bazhlekov, I. B. and Zapryanov, Z. D., "Interaction of a Deformable Bubble with a Rigid Wall at Moderate Reynolds Numbers," Journal of Fluid Mechanics, 219, pp. 241-271, 1990. Shyy,W., "A Study of Finite Difference Approximations to Steady-State, Convection-Dominated Flow Problems," Journal of Computational Physics, 57, pp. 415-438, 1985. Shyy, W., "A Note on Assessing Finite Difference Procedure for Large Peclet/Reynolds Number Flow Calculation," in: Boundarv and Interior Layers: Computational and Asy mptotic Methods , Miller, J. J. H. (ed.), 3, pp. 303-308, Boole Press, Dublin, 1984.

PAGE 207

190 Shyy, W. and Chen, M.-H., "Interaction of Thermocapillary and Natural Convection Flows During Solidification: Normal and Reduced Gravity Conditions." Journal of Crystal Growth, 108, pp. 247-261, 1991. Shyy, W. and Chen, M.-H., "Steady-State Natural Convection with Phase Change," International Journal of Heat and Mass Transfer, 33, pp. 2545-2563, 1990. Shyy, W., Pang, Y., Hunter, G. B., Wei, D. Y. and Chen, M.-H., "Modeling of Turbulent Transport and Solidification during Continuous Ingot Casting," International Journal of Heat and Mass Transfer, 35, pp. 1229-1245, 1992a. Shyy, W., Udaykumar, H. S. and Liang, S.-J., "An Interface Tracking Method Applied to Morphological Evolution During Phase Change, " AIAA 27 * Thermophysics Conference, AIAA 92-2902, 1992b. Silliman, W. J. and Scriven, L. E., "Separating Flow near a Static Contact Line: Slip at a Wall and Shape of a Free Surface, " Journal of Computational Physics, 34, pp. 287-313, 1980. Sparrow, E. M., Patankar, S. V. and Ramadhyani, S., "Analysis of Modeling in the Presence of Natural Convection in the Melt Region," ASME Journal of Heat Transfer, 99, pp. 520-526, 1977. Sullivan, J. M. Jr. and Lynch, D. R., "Non-Linear Simulation of Dendritic Solidification of an Undercooled Melt, " International Journal for Numerical Methods in Engineering, 25, pp. 415-444, 1988. Surek, T., Coriell, S. R. and Chalmers, B., "The Growth of Shaped Crystals from the Melt," Journal of Crystal Growth, 50, pp. 21-32, 1980. Swartz, J. C, Surek, T. and Chalmers, B., "The EFG Process Applied to the Growth of Silicon Ribbons," Journal of Electronic Materials, 4, pp. 255-279, 1975. Tatarchenko, V. A. and Brenner, E. A., "Crystallisation Stability during Capillarv Shaping, I and II, "Journal of Crystal Growth, 50, pp. 33-50, 1980. Thomas, P. D.,Ettouney, H. M. and Brown, R. A., "A thermal-Capillary Mechanism for a Growth Rate Limit in Edge-Defined Film-Fed Growth Rate Limit Silicon Sheets," Journal of Crystal Growth. 76, pp. 339-351, 1986. Thompson, J. F.,Warsi, Z. U. A. and Mastin, C. W.. Numerical Grid Generatinn Elsevier, New York, 1985.

PAGE 208

191 Thompson, M. E. and Szekely, J., "The Transient Behavior of Weldpools with a Deformed Free Surface," International Journal of Heat and Mass Transfer, 32, pp. 1007-1019, 1989. Tseng, A. A.,Zou, J., Wang, H. P. and Hoole, S. R., "Numerical Modeling of Macro and Micro Behaviors of Materials in Processing: A Review," Journal of Computational Physics, 102, pp. 1-17, 1992. Ungar, L.H. and Brown, R. A., "Cellular Interface Morphologies in Directional Solidification. The One-Sided Model," Physical Review B, 29, pp. 1367-1380, 1984. Viskanta, R., "Heat Transfer During Melting and Solidification of Metals," ASME Journal of Heat Transfer, 110, pp. 1205-1219,1988. Viskanta, R., "Natural Convection in Melting and Solidification," in: Natural Convection: Fundamentals and Applications . Kakac, S.,Aung, W. and Viskanta, R. (eds.), Hemisphere, Washington, D.C., pp. 845-877, 1985. Viskanta, R., "Phase change heat transfer," in: Solar Heat and Storage: Latent Heat Materials . Lane, G. A. (ed.), CRC Press, Boca Raton, FL, pp. 153-222, 1983. VoUer, V.R., "Implicit Finite-Difference Solutions of the Enthalpy Formulation of Stefan Problems," IMA Journal of Numerical Analysis, 5, pp. 201-214, 1985. Voller, V.R. and Cross, M., "Accurate Solutions of Moving Boundary Problems Using the Enthalpy Method," International Journal of Heat and Mass Transfer, 24, pp. 545-556, 1981. Voller, V. R., Swaminathan, C. R. and Thomas, B. G., "Fixed Grid Techniques for Phase Change Problems: A Review," International Journal for Numerical Methods in Engineering, 30, pp. 875-898, 1990. Wang, G.-X. and Matthys, E. F., "Numerical Modelling of Phase Change and Heat Transfer during Rapid Solidification Processes: Use of Control Volume Integrals with Element Subdivision," International Journal of Heat and Mass Transfer, 35, pp. 141-153, 1992. Warming, R. and Hyett, B., "The Modified Equation Approach to the Stability and Accuracy Analysis of Finite-Difference Methods," Journal of Computational Physics, 14, 159-179, 1974. Warming, R. F.,Kutler, P. and Lomax, H., "Secondand Third-Order Noncentral Difference Schemes for Nonlinear Hyperbolic Equations," AIAA Journal, 11, pp. 189-196, 1973.

PAGE 209

192 Wilson, D. G., Solomon, A. D. and Boggs, P. T., Moving Boundary Problems, Academic Press, New York, 1978. Wouters, P.. van Schaftingen, J. J., Crochet, M. J. and Geyling, F. T., "Numerical Simulation of the Horizontal Bridgman Growth. Part 3: Calculation of the Interface," International Journal for Numerical Methods in Fluids, 7, pp. 131-153, 1987. Young, G. W. and Davis, S. H., " A Plate Oscillating across a Liquid Surface: Effect of Contact Angle Hysteresis," Journal of Fluid Mechanics, 174, pp. 327-356, 1987. Zulehner, W., "Czochralski Growth of Silicon," Journal of Crystal Growth, 65, pp. 189-213, 1983.

PAGE 210

BIOGRAPHICAL SKETCH Shin-Jye Liang was bom on April 9, 1958, in Maoli, Taiwan, the Republic of China, and spent much of his early life there, enjoying the rivers for summertime swimming and hiking the mountain during winter. He graduated with a Bachelor of Science degree in hydraulic engineering from National Cheng Kung University in 1980 and earned his Master of Science degree in Civil Engineering from National Taiwan University in 1983. After two years of service in ROTC in Taiwan, he came to the United States and received his second Master of Science degree, in ocean engineering, from the University of Rhode Island in 1988. The author enrolled as a Ph.D. student in the Department of Aerospace Engineering, Mechanics, and Engineering Science at the University of Florida in 1989. He enjoys many sports and reading. 193

PAGE 211

I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree o^^Â»octor of Philosophy. Wei Shyy, Chairman Professor of Aerospace Engineering, Mechanics, and Engineering Science I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. Chen-Chi Hsu Professor of Aerospace Engineering, Mechanics, and Engineering Science I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. Ulrich H. Kurzweg Professor of Aerospace Engineering, Mechanics, and Engineering Science I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. Assistant Professor of Aerospace Engineering, Mechanics, and Engineering Science

PAGE 212

I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate, in scope and quality, as a dissertation for the degree of Doctor of Philosophy. Craig Saltiell Assistant Pre ssorxif Mechanical Engineering This dissertation was submitted to the Graduate Faculty of the College of Engineering and to the Graduate School and was accepted as partial fulfillment of the requirements for the degree of Doctor of Philosophy. May 1993 Winfred M. Phillips Dean, College of Engineering Madelyn M. Lockhart Dean, Graduate School

xml version 1.0 encoding UTF-8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd
INGEST IEID ER828U6WP_OKN3X8 INGEST_TIME 2015-04-22T21:51:45Z PACKAGE AA00029996_00001
AGREEMENT_INFO ACCOUNT UF PROJECT UFDC
FILES