
Citation 
 Permanent Link:
 http://ufdc.ufl.edu/AA00029996/00001
Material Information
 Title:
 A study of free and moving boundary problems involving thin crystal growth
 Creator:
 Liang, ShinJye, 1958
 Publication Date:
 1993
 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 )
nonfiction ( marcgt )
Notes
 Thesis:
 Thesis (Ph. D.)University of Florida, 1993.
 Bibliography:
 Includes bibliographical references (leaves 183192).
 General Note:
 Typescript.
 General Note:
 Vita.
 Statement of Responsibility:
 by ShinJye Liang.
Record Information
 Source Institution:
 University of Florida
 Holding Location:
 University of Florida
 Rights Management:
 The University of Florida George A. Smathers Libraries respect the intellectual property rights of others and do not claim any copyright interest in this item. This item may be protected by copyright but is made available here under a claim of fair use (17 U.S.C. Â§107) for nonprofit research and educational purposes. Users of this work have responsibility for determining copyright status prior to reusing, publishing or reproducing this item for purposes other than what is allowed by fair use or other copyright exemptions. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder. The Smathers Libraries would like to learn more about this item and invite individuals or organizations to contact the RDS coordinator (ufdissertations@uflib.ufl.edu) with any additional information they can provide.
 Resource Identifier:
 030010770 ( ALEPH )
29809577 ( OCLC )

Downloads 
This item has the following downloads:

Full Text 
A STUDY OF FREE AND MOVING BOUNDARY PROBLEMS
INVOLVING THIN CRYSTAL GROWTH
By
ShinJye 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, ChenChi 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
TABLE OF CONTENTS
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 EdgeDefined 1.3 Scope of the Present Work
1.4 Outline ...............
FilmFed Growth Technique . . . . . . . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . . . . . . .
2 INTERACTION OF TIME STEPPING AND CONVECTION
FOR UNSTEADY FLOW SIMULATION ...............
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 Quasiequilibrium 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 Quasistationary 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
21 Amplification Factors of Backward Euler Scheme for Linear
Burgers Equation ................................... 35
22 Amplification Factors of CrankNicolson Scheme for Linear
Burgers Equation ................................... .36
41 Assessment of time stepping scheme for phase change problem.
(a) Comparisons of the predicted interface location at r = 0.2 using backward Euler and CrankNicolson 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 CrankNicolson time stepping schemes for Eqs.
(4.3) and (4.4) with 11 and 21 grids ..................... .132
42 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
43 Definition and magnitude of key dimensionless parameters based on
values given in Table 41 (Case 1: St = 1, Case 2: St = 0.024)
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134
44 Two different scaling procedure and resulting nondimensional
energy equations and boundary conditions ................ 135
vi
LIST OF FIGURES
Figures Page
11 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
22 Comparison of solutions of unsteady linear Burgers equation
obtained using the backward Euler and CrankNicolson time
stepping schemes .................................... 38
23 Comparison of the backward Euler and CrankNicolson time
stepping schemes with four convection schemes for unsteady linear
Burgers equation using 81 nodal points, and P = 10. C = 1 . ... 39
24 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
25 The L2 error norm with respect to Ax of the nonlinear Burgers
equation with v = 5x103 and At/Ax = 0.5 ................ 41
26 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
27 The L2 error norm with respect to Ax of the nonlinear Burgers
equation with k = 5x105 and At/Ax = 0.5 ................ 43
vii
28 Effect of time stepping schemes on solution accuracy of the
nonlinear Burgers equation with v = 5x10 , Ax = 0.02, and
At = 0.05 ......................................... 44
29 Performance characteristics of the combined backward Euler and
firstorder upwind schemes with periodic boundary condition,
P = 103 , C = 0.5,and Ax = 0.05 ....................... 45
210 Performance characteristics of the combined backward Euler and
secondorder central difference schemes with periodic boundary
condition, P = 103 , C = 0.5,and Ax = 0.05 ............... .46
211 Performance characteristics of the combined backward Euler and
secondorder upwind schemes with periodic boundary condition,
P = 103 , C = 0.5,and Ax = 0.05 ....................... 47
212 Performance characteristics of the combined backward Euler and
secondorder QUICK schemes with periodic boundary condition,
P = 103 , C = 0.5,and Ax = 0.05 ....................... 48
213 Performance characteristics of the combined CrankNicolson and
firstorder upwind schemes with periodic boundary condition,
P = 103 , C = 0.5,and Ax = 0.05 ....................... 49
214 Performance characteristics of the combined CrankNicolson and
secondorder central difference schemes with periodic boundary
condition, P = 103, C = 0.5, and Ax = 0.05 ................ 50
215 Performance characteristics of the combined CrankNicolson and
secondorder upwind schemes with periodic boundary condition,
P = 103 , C = 0.5,and Ax = 0.05 ....................... 51
216 Performance characteristics of the combined CrankNicolson and
secondorder QUICK schemes with periodic boundary condition,
P = 103 , C = 0.5,and Ax = 0.05 ....................... 52
31 Schematic of meniscus corresponding to edgedefined filmfed
growth process ..................................... 79
32 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*.
33 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' .
34 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 .
35 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 .
36 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.)
37 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
38 Energy profiles for Bo = 102 , AP = I%, , = 700 and various
he/rb. Both upward and downward pulling cases are shown . ... 87
39 Schematic presentation of hysteresis model ................ 88
310 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.
311 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.
312 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.
313 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.
314 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.
315 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
41 Exact solution for melting of a solid confined to a semiinfinite
region ........................................... 137
42 Comparison of interface trajectories predicted by the coupled
(full equations) and decoupled (quasistationary equations)
approaches with three values of the Stefan numbers ......... .138
43 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
44 Schematic configurations and boundary conditions of the present
EFG process ...................................... 140
45 Comparison of linear meniscus profile with the LaplaceYoung
solution with RI = 1, H, = 0.188, R, = 0.763, and various values of Bond number. (a) Ap = 5 % and (b) Ap = 0 % ........... 141
46 The grid system and isothermal contours of the initial conditions for
both Stefan numbers ................................ 142
47 The grid system and isothermal contours of the steadystate solution
for Case 1: St = I .................................. 143
48 The grid system and isothermal contours of the steadystate 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 = O(St) and choice 2:
VN = (1) ....................................... 145
51 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
52 Time histories of Hc and R, subject to a threeharmonic
perturbation of 0. with Q, = 500 Ar, '2 = 292 Ar, and
03 = 108 Ar (St = 1.0) .............................. . 164
53 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
54 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
55 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
56 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
57 Time histories He and Rc subject to a threeharmonic perturbation
of U, with Q, = 500 AT, (2, = 292 AT and (2 = 108 Ar (St = 1.059
58 Time histories of sawtoothed 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).
59 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
510 Characteristics of interface movement with respect to different
forcing periods of U, (St = 0.024) ...................... 172
511 A typical Bode diagram depicting sensitivity of the fibre diameter
variation with respect to different forcing frequencies of U,
(St = 0.024) ...................................... 173
512 Model of dynamic contact condition at the trijunction point
(St = 0.024) ...................................... 174
513 Time histories of H. and R, subject to a single harmonic
perturbation of U, with period = L00Ar (St = 0.024) ........ 175
514 Time histories of dH/dr and dRlIdr subject to a single
harmonic perturbation of U, with period = LOOAT (St = 0.024) 176
515 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 AB/ 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 mth component thermal conductivity nondimensional 1 wavelength of the mth 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
AM AdamsBashforth/AdamsMoulton predictorcorrector scheme
Backward Euler time stepping scheme
Biot number
Bond number
Courant number
CrankNicolson 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
L2 error norm
process EdgeDefined Filmfed 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
VI
[cm']
[W/cmK]
[cm] [cm] [J/g] [cm]
xiii
[cm]
[J/gK]
[erg] [erg] [erg] [erg]
[cm]
W/cm2] [cm/s2]
[cm]
cm2K]
[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. (444) nondimensional radius
radius [cm]
Rayleigh number
Radiation number
radius of die [cm]
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. (410) dynamic viscosity [g/
kinematic viscosity [c
transformed coordinates density [g
StefanBoltzman constant [5.670x 1012 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 steadystate apparent contact angle period of external forcing of the i* harmonics
icm] m2/s]
/cm3] 2K4]
[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
n1
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
ShinJye 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 Edgedefined Filmfed 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 quasidynamic 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 Stefantype (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 edgedefined filmfed growth (EFG) process for producing thin sapphire crystals.
1.2 Description of Edgedefined Filmfed Growth Technique
Sapphire fibers for use in both optical sensors and structural composites have been produced using the edgedefined filmfed 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 11 (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 watercooled, environmentally regulated chamber. The melt is supplied to each die tip from a capillaryfed 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 wellcontrolled 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 shapedefining dies can be enclosed within a watercooled, environmentally regulated chamber, as illustrated in Figure 11. 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 12. 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 timedependent 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 timedependent 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 solidliquidgas
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 liquidgassolid 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. Bodyfitted 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 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. Quasiequilibrium 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 singlevalued 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 11 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
convection & radiation
Figure 12
pulling speed
Solid
(diffusion)
interface
trijunction (latent heat) point
meniscus Liquid (YoungLaplace 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 firstorder upwind scheme used to be a popular choice due to its superior numerical stability as well as the wigglesfree characteristics. However, the resulting accuracy yielded by this scheme is often found to be inadequate in view of its firstorder 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 convectiondominated 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 firstorder (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 firstorder backward Euler, the secondorder CrankNicolson, and the thirdorder AdamsBashforth/AdamsMoulton (AB/AM) predictorcorrector schemes, are considered. In conjunction, four convection schemes, firstorder upwind, secondorder central differencing, secondorder 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 43.
2.2 Model Problem Analysis
The model problem considered is the unsteady onedimensional 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 secondorder 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) firstorder 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) secondorder central differencing
a(u(0) _ (uk)1,  (u)i 1 + T (2.5)
ax 2 Ax
where T, is O[(Ax)2]
(c) secondorder upwind
1 [ 3(uO), 4(uo),_ +(uo)i2 ] + 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 firstorder backward Euler method for the time derivative is used, the discretized form of Eq. (2. 1) can be expressed as follows:
(a) firstorder 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) secondorder central differencing
P P P P 0 29
(1+f) 4,_1  (2+) 4, + (1 ) 4,, =   4 (2.9)
2 C' 2 C'
(c) secondorder 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 CrankNicolson scheme is a secondorder 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) firstorder upwind
(1+P) 4  (2+P+21) P + 4k., = C (2.12)
 (1 +P) q5 + (2+P2) 4
C
(b) secondorder central differencing
(1+P) 4,_  (2+2)4 +(1 ) 4,) =
2 C 2 (2.13)
 (1+) 0_ + (22.) P v  (l.E) 0.
2 C ' 2 '+1
21
(c) secondorder upwind
 Oi2 +(I+2P) Oi  (2+3P+2) O, + 4,=
2 2 C (2.14)
P (i + (2+3P2P) 4  .
2 2 C
(d) QUICK
P2 + (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 + P2 ) 00  0 P) i+.
8 8 8 C 8
The thirdorder predictorcorrector scheme, based on the AdamsBashforth scheme as the predictor step and the AdamsMoulton 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(P4) =J(x,k) (2.16)
at ax ax ax
then the present predictorcorrector algorithm results AdamsBashforth method :
0(00 = O(n) + A( 3f(")  f("1 ) (2.17)
22
AdamsMoulton method :
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 AB/AM 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 CrankNicolson 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 firstorder and secondorder upwind solutions are wigglesfree while neither secondorder central differencing nor QUICK is wigglesfree (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 AB/AM 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 21 and 22 for the two time stepping schemes at t = 0.3,0.8, and 2. Using the backward Euler scheme, as shown in Figure 21, both the firstorder and secondorder upwind solutions are wigglesfree, but the firstorder 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 21 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.
At t = 2, as shown in Figure 21 (c), the solutions have approached steady state. Both the first and secondorder upwind schemes converge to virtually the same solution. Even though the firstorder upwind scheme works nearly as well as the secondorder upwind convection scheme for the onedimensional steady state equation, when combined with the backward Euler time stepping scheme, it is less satisfactory for the unsteady equation in view of the stronger numerical smearing effects depicted in Figures 21 (a) and (b). For the other two convection schemes, while the secondorder central difference scheme does not perform well for either steady or unsteady problems with strong convection effect, QUICK is more satisfactory for the unsteady equation, yielding wiggles mainly when the steady state is reached. This is because at steady state the solution variation is confined in a thin boundary layer close to the right hand boundary, resulting in a higher gradient there. Overall, the secondorder upwind scheme performs best if both steady and unsteady solutions are needed.
Figure 22 compares the solutions obtained by using the backward Euler and CrankNicolson schemes; they are noticeably different. For the firstorder upwind scheme, with C = 1, the initial profile is better preserved and the gradient is less smeared out with CrankNicolson. In conjunction with the firstorder upwind
25
scheme, CrankNicolson, being less dissipative but yielding no wiggles, performs better than backward Euler. Similarly, the performance of a secondorder time stepping scheme such as CrankNicolson is less satisfactory when combined with other secondorder convection schemes due to the fact that both time stepping and convection schemes are dispersive. As illustrated in Figures 22 (b) (d), nonphysical solutions are observed for all three secondorder convection schemes. For both the central differencing and QUICK schemes, the oscillations yielded by CrankNicolson are substantially larger than those by backward Euler. For the secondorder upwind scheme, while the CrankNicolson 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 21 and 22 are generic and not affected by the variations in the number of grid points. In order to demonstrate this fact, Figure 23 compares the performances between the CrankNicolson 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 23 can be observed in Figure 22.
2.3.2Nonlinear Burgers Equation
For the nonlinear Burgers equation, two different values of viscosity, P =
5x103 and 5x105, are used to explore the effect of solution gradients on the
26
performance of each scheme. Three time stepping schemes, including the firstorder backward Euler scheme, the secondorder CrankNicolson scheme, and the thirdorder AdamsBashforth/AdamsMoulton predictorcorrector scheme, are studied along with the firstorder and secondorder 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 24 compares the solutions yielded by the three time stepping schemes in conjunction with (a) the firstorder upwind convection scheme and (b) the secondorder upwind convection scheme. All the computations shown in Figure 24 are based on 51 nodes (Ax = 0.02), and At = 0.01, resulting in P = 4u and C = 0.5u. Figure 24 shows that with firstorder upwind, the numerical solutions yielded by backward Euler and CrankNicolson are fairly comparable while noticeable improvement can be obtained by using the thirdorder AB/AM scheme. Figure 24 indicates that a firstorder scheme in space and a thirdorder scheme in time can actually form a good combination.
Regarding the use of the secondorder upwind convection scheme, shown in Figure 24 (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 25 depicts an L2 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 L2 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 25. 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 firstorder upwind convection scheme is proportional to Ax, Figure 25 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 25 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 firstorder upwind convection alone. As to the secondorder upwind convection scheme, Figure 25 (b) shows that the solutions using both CrankNicolson and AB/AM 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 secondorder upwind convection scheme, they are essentially proportional to (Ax)2. Overall, Figures 24 and 25 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 26 compares the numerical solutions corresponding to the different combinations of the convection and time stepping schemes. Comparing Figures 24 and 26, 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 26 (a) shows that, even with the use of the firstorder upwind convection scheme, which is viewed as usually too dissipative, nonphysical wiggles develop when the AB/AM time stepping scheme is employed. With either backward Euler or CrankNicolson, 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 26 (b), in which the numerical wiggles appear in all solutions.
The L2 error norms of the solutions of v = 5x105 are summarized in Figure 27 where the Courant number again remains unchanged as Ax is varied. Unlike those curves shown in Figure 25 for v = 5x103, the results in Figure 27 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 AB/AM 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 higherorder 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 28, where both v, 5x105, and Ax, 0.02, remain the same as those used in Figure 26, 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 higherorder time stepping schemes exhibit stronger dispersive characteristics. Hence, as demonstrated in Figure 28 the firstorder upwind convection scheme, being generally more diffusive, can be used along with the CrankNicolson 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
secondorder 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 = fl,
km = 27r/L m is the wave number of the mth component,
and Lm is the wavelength of the mth 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 CrankNicolson 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 29 to 216 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 29 to 212; the corresponding results of the CrankNicolson time stepping scheme are presented in Figures 213 to 216.
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 29 to 216. First, CrankNicolson 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 29 to 216 demonstrate that for all three upwind type of convection schemes, CrankNicolson 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 secondorder central difference scheme are very close to the exact values; when combined with CrankNicolson, 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 secondorder upwind scheme is actually more dissipative than the firstorder upwind scheme at short wavelengths as indicated by part (iii) of Figures 29 and 211. With both CrankNicolson and backward Euler, the secondorder upwind scheme is less dissipative than the firstorder 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 211 and 215, it becomes clear that it is the phase errors in the long wavelength that are mostly responsible for the inaccuracies created by the combined CrankNicolson/secondorder upwind schemes. This can explain why its solutions exhibit a single long wave overshoot instead of 2Ax short wave oscillations. For the firstorder 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/secondorder upwind or CrankNicolson/firstorder 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, CrankNicolson 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 secondorder upwind scheme is less dissipative than the firstorder upwind scheme mainly for the longer waves. Consequently, for the secondorder upwind scheme, backward Euler yields more accurate phase speed than CrankNicolson; while for the firstorder upwind scheme, it is CrankNicolson that is more satisfactory in phase speed. Regarding the secondorder 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 secondorder 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 firstorder upwind for convection and CrankNicolson for time, or a combination of secondorder upwind for convection and backward Euler for time performs better. A consistent secondorder 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 21 Amplification Factors of Backward Euler Scheme
for Linear Burgers Equation
scheme amplification factor
firstorder p
upwind scheme
G = C
[(2+P+) + (2+P)cos ]  i(Psin )
C
secondorder p
central difference
[(2+) + 2cosp I  i(Psinp)
C
secondorder 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 22 Amplification Factors of CrankNicolson Scheme for Linear Burgers Equation
scheme amplification factor
firstorder
upwind scheme [(2 +P2)  (2+P)cosp] + i(Psins)
G5 C
(2 +P+2) + (2+P)cos p I  i(Psinp)
C
secondorder P
central difference [(22)  2cosp] + i(Psin p)
G6 C
[(2+) + 2cosp I  i(Psinp)
C
secondorder 3 P P P
upwind scheme [(2+P2) + 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+P2)+ 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 secondorder central differencing scheme
3 secondorder 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 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.
2
0 0.2 0.4
C
0
b
V
1.2
1.0 0.8 0.6
0.4 0.2
0
(a) firstorder 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
secondorder central differencing scheme d. =0.8 1=03
1.2
(b) secondorder upwind scheme (d) QUICK scheme
Figure 22 Comparison of solutions of unsteady linear Burgers equation obtained using
the backward Euler and CrankNicolson 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) secondorder upwind scheme
, =0.8
, /t = .3
I. /
a
(a) firstorder upwind scheme
0
0.2 0.4 0.6 0.8
x
1.0
1.2
(d) QUICK scheme
Figure 23 Comparison of the backward Euler and CrankNicolson time stepping
with four convection schemes for unsteady linear Burgers equation
nodal points, and P = 10, C = 1.
schemes using 81
108
izO
c ~
(c) secondorder 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) firstorder upwind convection scheme
(b) secondorder upwind convection scheme
Figure 24 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 AB/AM
I i iI
a
0
0.8 1.0 1.2
10.2
I 4
10i 102 10
t 2
2 3
1 Backward Euler 2 Crank Nicolson
3 AB/AM
= 3.5
2
,ag3
104 [
105
101 10
Ax
t=2
2
3
t =3.s
3
1 102 10
b Ax
(a) firstorder upwind convection scheme
(b) secondorder upwind convection scheme
Figure 25 The L2 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 A8/AM Exact
solution
1 2
0.2 0.4 0.6 0.8 1.0
X
(a) firstorder upwind convection scheme
(b) secondorder upwind convection scheme
Figure 26 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
101 10.2
b
(a) firstorder upwind convection scheme
Figure 27
(b) secondorder upwind convection scheme
The L2 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 CrankNicolson
3 AB/AM
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) firstorder upwind convection scheme (b) secondorder upwind convection scheme
Figure 28 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
3AB/AM
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 29 Performance characteristics of the combined backward Euler and firstorder
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 210
Performance characteristics central difference schemes C = 0.5,and Ax = 0.05.
of the combined backward Euler and secondorder 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 211
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 secondorder 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 212 Performance characteristics of the combined backward Euler and secondorder
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 213 Performance characteristics of the combined CrankNicolson and firstorder
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 214
Performance characteristics of the combined CrankNicolson 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 215
Performance characteristics order upwind schemes with and Ax = 0.05.
of the combined CrankNicolson 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 216 Performance characteristics of the combined CrankNicolson and secondorder
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 crystalpulling 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 liquidgas 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 viceversa, 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 LaplaceYoung equation for the equilibrium meniscus profiles. The boundary conditions are chosen to correspond to that of the edgedefined filmfed 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 31, 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 LaplaceYoung equation for the liquidgas 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 coordinate defining the surface of the meniscus, yi is the surface tension of the liquidgas 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 4x103 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 nondimensionalize all physical quantities as follows, *'s representing nondimensional 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 nondimensionalized profile equation, after having dropped all *'s is,
____ 1
3 ZAp (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 nondimensionalized 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 solidgas, solidliquid and liquidgas interface. , is the static contact angle measured with respect to the liquidsolid 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 threephase 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 liquidgas 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 LaplaceYoung 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 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 liquidvapor 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 nonunique 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 = (2yys) 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 LaplaceYoung equation, is nonlinear and has to be solved numerically. The integration of the nonlinear ordinary differential Eq. (3.4) is performed using the fourthorder Adams predictorcorrector method, using the RungeKutta 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 Edgedefined Filmfed 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 32. 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 LaplaceYoung equation. Also shown in Figure 32 are sample meniscus profiles, corresponding to 0(, of 60* and 900, respectively.
As shown in Figure 32, 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 nonextremum 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 nonextremum may all be observed at the given contact angle 0.. 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.
65
As discussed in connection with the nondimensional 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 LaplaceYoung equation fail to exist.
Cases 2 and 3 (Figs. 33 and 34 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.5x106,Ap=5%, he/r. =2. For this low Bond number case, from Figure 33, 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 LaplaceYoung equation) obtained.
The contrasting situation to Case 2 emerges in Case 3 (Figure 34). 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 33 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 nondimensionalization. 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 zaxis.
67
Another aspect to note in Figs. 33 and 34 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 35). 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 35
(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 LaplaceYoung 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 36 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 LaplaceYoung 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 36.
69
With Bo = 102, when hC/rb 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 0, at the top dictates the meniscus shape. Figure 36 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 37 depicts the solution characteristics for the case of Bo = 102, 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 37 are presented in Figure 38; 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.3Ouasiequilibrium 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 LaplaceYoung 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 39. 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 quasiequilibrium 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 04, is the deviation in contact angle from the value that leads to fibre growth with nonuniform 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 nonlinearity 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 310, 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 310 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 311, 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. 311 (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 311 can be interpreted as a stickslip motion (Jansons 1986).
When we reach a higher Bond number, say 2.5x103 as in Figure 312, 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 310, 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 310, 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 310, to convex. For a higher meniscus, hc* = 1.5, for the same 40 and Ap shown in Figure 313, the radius variation is also not very substantial in magnitude and the waveform is not identical to hc(t), indicating the nonlinearity that is introduced by 4(t). When hysteresis is present, and for the case of fairly large hysteresis, and tall meniscus, Figure 314 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 314 (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 quasiequilibrium 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 quasiequilibrium situation, and the full fluiddynamic
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 fourthorder 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 nonextrema 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 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 nonequilibrium 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. Quasiequilibrium 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 31 Schematic of meniscus corresponding to edgedefined 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 32
81
2 1.5 1 0.5 0 0.5 1 1.5 2
f/
A0
4.5
4.
3.5
3,
Z5.
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
26
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) nonexI2MUM (minimum
0M 40 60 80 100 120 140
Figure 32  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) nonexurmum (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 33 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) nonCxUemUm
(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 34
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 nonexmum
2

Full Text 
PAGE 1
A STUDY OF FREE AND MOVING BOUNDARY PROBLEMS INVOLVING THIN CRYSTAL GROWTH By ShinJye 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, ChenChi 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 EdgeDefined FilmFed 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 Quasiequilibrium 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 Quasistationary 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 21 Amplification Factors of Backward Euler Scheme for Linear Burgers Equation 35 22 Amplification Factors of CrankNicolson Scheme for Linear Burgers Equation 36 41 Assessment of time stepping scheme for phase change problem. (a) Comparisons of the predicted interface location at r = 0.2 using backward Euler and CrankNicolson 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 CrankNicolson time stepping schemes for Eqs. (4.3) and (4.4) with 11 and 21 grids 132 42 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 43 Definition and magnitude of key dimensionless parameters based on values given in Table 41 (Case 1: St = 1, Case 2: St = 0.024) 134 44 Two different scaling procedure and resulting nondimensional energy equations and boundary conditions 135 vi
PAGE 7
LIST OF FIGURES Fi gures Page 11 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 22 Comparison of solutions of unsteady linear Burgers equation obtained using the backward Euler and CrankNicolson time stepping schemes 38 23 Comparison of the backward Euler and CrankNicolson time stepping schemes with four convection schemes for unsteady linear Burgers equation using 81 nodal points, and P = 10, C = 1 .... 39 24 Effect of time stepping schemes on solution accuracy of the nonlinear Burgers equation with y = 5x10'^ , Ax = 0.02, and At = 0.01 40 25 The L2 error norm with respect to Ax of the nonlinear Burgers equation with p = 5x10"^ and At/ Ax =0.5 41 26 Effect of time stepping schemes on solution accuracy of the nonlinear Burgers equation with u = 5x10'^ , Ax = 0.02. and At = 0.01 42 27 The L2 error norm with respect to Ax of the nonlinear Burgers equation with i> = 5x10'^ and At/ Ax =0.5 43 vii
PAGE 8
28 Effect of time stepping schemes on solution accuracy of the nonlinear Burgers equation with v = 5x10 ' , Ax = 0.02, and At = 0.05 44 29 Performance characteristics of the combined backward Euler and firstorder upwind schemes with periodic boundary condition, p = 10* , C = 0.5, and Ax = 0.05 45 210 Performance characteristics of the combined backward Euler and secondorder central difference schemes with periodic boundary condition, P = 10^ , C = 0.5, and Ax = 0.05 46 211 Performance characteristics of the combined backward Euler and secondorder upwind schemes with periodic boundary condition, P = 10' , C = 0.5,and Ax = 0.05 47 212 Performance characteristics of the combined backward Euler and secondorder QUICK schemes with periodic boundary condition, P = 10* , C = 0.5, and Ax = 0.05 48 213 Performance characteristics of the combined CrankNicolson and firstorder upwind schemes with periodic boundary condition, P = 10\ C = 0.5, and Ax = 0.05 49 214 Performance characteristics of the combined CrankNicolson and secondorder central difference schemes with periodic boundary condition, P = 10\ C = 0.5, and Ax = 0.05 50 215 Performance characteristics of the combined CrankNicolson and secondorder upwind schemes with periodic boundary condition, P = 10\ C = 0.5, and Ax = 0.05 51 216 Performance characteristics of the combined CrankNicolson and secondorder QUICK schemes with periodic boundary condition, P = 10\ C = 0.5, and Ax = 0.05 52 31 Schematic of meniscus corresponding to edgedefined filmfed growth process 79 32 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Â° . 33 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Â° . 34 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Â° . 36 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.) 37 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
38 Energy profiles for Bo = 10'^ , AP = 1%, 0o = 70Â° and various hc/rfc. Both upward and downward pulling cases are shown .... 87 39 Schematic presentation of hysteresis model 88 310 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. 311 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. 312 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. 313 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. 314 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. 315 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
41 Exact solution for melting of a solid confined to a semiinfinite region 137 42 Comparison of interface trajectories predicted by the coupled (full equations) and decoupled (quasistationary equations) approaches with three values of the Stefan numbers 138 43 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 44 Schematic configurations and boundary conditions of the present EFG process 140 45 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 46 The grid system and isothermal contours of the initial conditions for both Stefan numbers 142 47 The grid system and isothermal contours of the steadystate solution for Case 1: St = 1 143 48 The grid system and isothermal contours of the steadystate 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 52 Time histories of He and subject to a threeharmonic perturbation of with fi, = 500 At, U2 = 292 At, and ^3 = 108 At (St = 1.0) 164 53 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
54 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 55 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 56 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 57 Time histories H, and R^ subject to a threeharmonic perturbation of Up with fi, = 500 At, = 292 At and fij = 108 At (St = l.Q)69 58 Time histories of sawtoothed 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). 59 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 510 Characteristics of interface movement with respect to different forcing periods of Up (St = 0.024) 172 511 A typical Bode diagram depicting sensitivity of the fibre diameter variation with respect to different forcing frequencies of Up (St = 0.024) 173 512 Model of dynamic contact condition at the trijunction point (St = 0.024) 174 513 Time histories of and R^ subject to a single harmonic perturbation of Up with period = IOOAt (St = 0.024) 175 514 Time histories of dH^/dr and dR^/dr subject to a single harmonic perturbation of Up with period = IOOAt (St = 0.024) 176 515 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] AB/AM AdamsBashforth/AdamsMoulton predictorcorrector scheme B. E, Backward Euler time stepping scheme Hi Biot number t , Bo Bond number C Courant number C. N. CrankNicolson time stepping scheme Cp specific heat at constant pressure [J/gK] 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 L2 error norm EFG process EdgeDefined Filmfed 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 Vl ij integers denoting grid point in space coordinate Ka wavenumber of the mth component [cm '] k thermal conductivity [W/cmK] L nondimensional 1 wavelength of the mth 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. (444) 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. (410) fi dynamic viscosity [g/scm] u kinematic viscosity [cm^/s] transformed coordinates p density [g/cm^] a StefanBoltzman constant [5.670xlO'^W/cmK'*] 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 steadystate [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 n1 previous time step n current time step nf1 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 ShinJye 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 Edgedefined Filmfed 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 quasidynamic 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 Stefantype (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 edgedefmed filmfed growth (EFG) process for producing thin sapphire crystals. 1.2 Description of Edgedefined Filmfed Growth Technique Sapphire fibers for use in both optical sensors and structural composites have been produced using the edgedefined filmfed 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 11 (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 watercooled, environmentally regulated chamber. The melt is supplied to each die tip from a capillaryfed 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 wellcontrolled 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 shapedefining dies can be enclosed within a watercooled, environmentally regulated chamber, as illustrated in Figure 11. 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 12. 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 timedependent 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 timedependent 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 solidliquidgas 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 liquidgassolid 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. Bodyfitted 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. Quasiequilibrium 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 singlevalued 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 (YoungLaplace Eq) Figure 12 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 firstorder upwind scheme used to be a popular choice due to its superior numerical stability as well as the wigglesfree characteristics. However, the resulting accuracy yielded by this scheme is often found to be inadequate in view of its firstorder 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 convectiondominated 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 firstorder (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 firstorder backward Euler, the secondorder CrankNicolson. and the thirdorder AdamsBashforth/AdamsMoulton (AB/AM) predictorcorrector schemes, are considered. In conjunction, four convection schemes, firstorder upwind, secondorder central differencing, secondorder 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 43. 2.2 Model Problem Analysis The model problem considered is the unsteady onedimensional 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 secondorder 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) firstorder 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) secondorder central differencing d(u({>) I ^ iu\., ^ dx 2Ax where T, is 0[(Ax)^] (c) secondorder 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 firstorder backward Euler method for the time derivative is used, the discretized form of Eq. (2.1) can be expressed as follows: (a) firstorder upwind where P = u Ax/p = cell Peclet number, C = u At/ Ax = Courant number and the superscript "o" represents the known value. (b) secondorder central differencing (c) secondorder upwind
PAGE 37
20 (d) QUICK (2.11) The CrankNicolson scheme is a secondorder 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) firstorder upwind (2.12) (b) secondorder central differencing (Ul) 0,.. (22Z) 0,. .(1^) 0,.^, = (2.13)
PAGE 38
(c) secondorder upwind 21 (2.14) (d) QUICK (2.15) The thirdorder predictorcorrector scheme, based on the AdamsBashforth scheme as the predictor step and the AdamsMoulton 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 predictorcorrector algorithm results : AdamsBashforth method : (2.17)
PAGE 39
AdamsMoulton 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 AB/AM 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 CrankNicolson 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 firstorder and secondorder upwind solutions are wigglesfree while neither secondorder central differencing nor QUICK is wigglesfree (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 AB/AM 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 21 and 22 for the two time stepping schemes at t = 0.3,0.8, and 2. Using the backward Euler scheme, as shown in Figure 21, both the firstorder and secondorder upwind solutions are wigglesfree, but the firstorder 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 21 shows that this scheme also produces nonphysical oscillations behind the sharp front for the unsteady
PAGE 41
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. At t = 2, as shown in Figure 21 (c), the solutions have approached steady state. Both the firstand secondorder upwind schemes converge to virtually the same solution. Even though the firstorder upwind scheme works nearly as well as the secondorder upwind convection scheme for the onedimensional steady state equation, when combined with the backward Euler time stepping scheme, it is less satisfactory for the unsteady equation in view of the stronger numerical smearing effects depicted in Figures 21 (a) and (b). For the other two convection schemes, while the secondorder central difference scheme does not perform well for either steady or unsteady problems with strong convection effect, QUICK is more satisfactory for the unsteady equation, yielding wiggles mainly when the steady state is reached. This is because at steady state the solution variation is confined in a thin boundary layer close to the right hand boundary, resulting in a higher gradient there. Overall, the secondorder upwind scheme performs best if both steady and unsteady solutions are needed. Figure 22 compares the solutions obtained by using the backward Euler and CrankNicolson schemes; they are noticeably different. For the firstorder upwind scheme, with C = 1, the initial profile is better preserved and the gradient is less smeared out with CrankNicolson. In conjunction with the firstorder upwind
PAGE 42
25 scheme, CrankNicolson, being less dissipative but yielding no wiggles, performs better than backward Euler. Similarly, the performance of a secondorder time stepping scheme such as CrankNicolson is less satisfactory when combined with other secondorder convection schemes due to the fact that both time stepping and convection schemes are dispersive. As illustrated in Figures 22 (b) (d), nonphysical solutions are observed for all three secondorder convection schemes. For both the central differencing and QUICK schemes, the oscillations yielded by CrankNicolson are substantially larger than those by backward Euler. For the secondorder upwind scheme, while the CrankNicolson 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 21 and 22 are generic and not affected by the variations in the number of grid points. In order to demonstrate this fact, Figure 23 compares the performances between the CrankNicolson 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 23 can be observed in Figure 22. 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 firstorder backward Euler scheme, the secondorder CrankNicolson scheme, and the thirdorder AdamsBashforth/AdamsMoulton predictorcorrector scheme, are studied along with the firstorder and secondorder 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 24 compares the solutions yielded by the three time stepping schemes in conjunction with (a) the firstorder upwind convection scheme and (b) the secondorder upwind convection scheme. All the computations shown in Figure 24 are based on 51 nodes (Ax = 0.02), and At = 0.01, resulting in P = 4u and C = 0.5u. Figure 24 shows that with firstorder upwind, the numerical solutions yielded by backward Euler and CrankNicolson are fairly comparable while noticeable improvement can be obtained by using the thirdorder AB/AM scheme. Figure 24 indicates that a firstorder scheme in space and a thirdorder scheme in time can actually form a good combination. Regarding the use of the secondorder upwind convection scheme, shown in Figure 24 (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 25 depicts an L2 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 L2 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 26 compares the numerical solutions corresponding to the different combinations of the convection and time stepping schemes. Comparing Figures 24 and 26, 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 26 (a) shows that, even with the use of the firstorder upwind convection scheme, which is viewed as usually too dissipative, nonphysical wiggles develop when the AB/AM time stepping scheme is employed. With either backward Euler or CrankNicolson, 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 26 (b), in which the numerical wiggles appear in all solutions. The L2 error norms of the solutions of u = 5x10'* are summarized in Figure 27 where the Courant number again remains unchanged as Ax is varied. Unlike those curves shown in Figure 25 for v = 5x10'^, the results in Figure 27 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 AB/AM 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 higherorder 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 28, where both 5xl0^and Ax, 0.02, remain the same as those used in Figure 26, 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 higherorder time stepping schemes exhibit stronger dispersive characteristics. Hence, as demonstrated in Figure 28 the firstorder upwind convection scheme, being generally more diffusive, can be used along with the CrankNicolson 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 secondorder 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 = V1, kÂ„ = 2ir/LÂ„ is the wave number of the mth component, and LÂ„ is the wavelength of the mth component. Based on the elementary solution, the amplification factor of the exact solution can be defined
PAGE 48
31 G = ^(^^^) = e'^^'e'^"" (221) 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 CrankNicolson 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 29 to 216 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 29 to 212; the corresponding results of the CrankNicolson time stepping scheme are presented in Figures 213 to 216.
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 29 to 216. First, CrankNicolson 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 29 to 216 demonstrate that for all three upwind type of convection schemes, CrankNicolson 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 secondorder central difference scheme are very close to the exact values; when combined with CrankNicolson, 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 secondorder upwind scheme is actually more dissipative than the firstorder upwind scheme at short wavelengths as indicated by part (iii) of Figures 29 and 211. With both CrankNicolson and backward Euler, the secondorder upwind scheme is less dissipative than the firstorder 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 211 and 215, it becomes clear that it is the phase errors in the long wavelength that are mostly responsible for the inaccuracies created by the combined CrankNicolson/secondorder upwind schemes. This can explain why its solutions exhibit a single long wave overshoot instead of 2Ax short wave oscillations. For the firstorder 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/ secondorder upwind or CrankNicolson/firstorder 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, CrankNicolson 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 secondorder upwind scheme is less dissipative than the firstorder upwind scheme mainly for the longer waves. Consequently, for the secondorder upwind scheme, backward Euler yields more accurate phase speed than CrankNicolson; while for the firstorder upwind scheme, it is CrankNicolson that is more satisfactory in phase speed. Regarding the secondorder 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 secondorder 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 firstorder upwind for convection and CrankNicolson for time, or a combination of secondorder upwind for convection and backward Euler for time performs better. A consistent secondorder 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 + aTu 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 CrankNicolson 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 wi 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 . Â— ^ CO il Â« o u ^
PAGE 70
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 infmite extent. On the other hand, in the crystalpulling 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
PAGE 71
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 (f>^ 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 4>^ is greater than that formed by receding liquid, say (f)^. Thus, the apparent angle assumed under static condition by the liquidgas surface at the trij unction is either or <^r . 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 to <^R or viceversa, depending on the direction of impending motion.
PAGE 72
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 f^ 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 LaplaceYoung equation for the equilibrium meniscus profiles. The boundary conditions are chosen to correspond to that of the edgedefined filmfed 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.
PAGE 73
56 3.2 Configuration. Assumptions and Formulation The configuration chosen here is shown in Figure 31, where hÂ„ re, and rt 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 LaplaceYoung equation for the liquidgas interface is applied, which in dimensional form is (Surek et al. 1980, Pitts 1976, and Huh and Scriven 1960) = Ap5z Ap where f is the radial coordinate defining the surface of the meniscus, is the surface tension of the liquidgas interface, Ap is the difference in density between liquid and gas phases, and Ap is the difference in applied pressure between the liquid, pi , and gas phases, Pg , i.e.,Ap = Pi 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
PAGE 74
57 a= JJl. (3.2) where, for example, for sapphire (AI2O3), 7,g = 0.69 N/m, Ap = 3.4x10? kg/m^ , and hence the capillary length is about 4x10" m. the bond number, which relates radius of base or die to capillary length, and controls relative importance gravity forces is given by , bo="^" (33) a' we specify aspect ratio, hjty, bo, each case. meniscus height calculated from ratio r^. nondimensionalize all physical quantities as follows,="^'s" representing nondimensional values, z*="z/a" ,r="f/a" ap*="Ap/(Apga)," ap '="1." again, using sapphire an example, reference pressure p,="Ap" g a about 130 n m^ thus, nondimensionalized profile equation, after having dropped *'s is, fzz 1 t r="^^p" (3.4) this equation seemingly invariant with respect since number does not explicitly appear eq. (3.4). misleading impression, influence profiles via term z, term. for h, rb, effect
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 solidgas, solidliquid and liquidgas interface. is the static contact angle measured with respect to the liquidsolid 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 threephase 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 liquidgas 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 liquidvapor 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 nonunique 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 (38) 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 fourthorder Adams predictorcorrector method, using the RungeKutta 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 Edgedefined Filmfed 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 32. 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 32 are sample meniscus profiles, corresponding to (f>o of 60Â° and 90Â°, respectively. As shown in Figure 32, 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 nonextremum 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 nonextremum 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 nondimensional 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. 33 and 34 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 33, 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 34). 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 33 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 nondimensionalization. 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 zaxis.
PAGE 84
67 Another aspect to note in Figs. 33 and 34 is that in the E vs. ^ curves both maxima and minima are obtained at specified Â„ = 130Â°, Figures 35 (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 36 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 36.
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 36 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 37 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 37 are presented in Figure 38; 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 Ouasiequilibrium 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 39. 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,iy,g) cannot be replaced by an expression such as 7rre^7,gCos<^e. for otherwise (7,i7,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)^ (313) where (f>o is the deviation in contact angle from the value that leads to fibre growth with nonuniform 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 nonlinearity 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 310 for the imposed amplitude AhÂ„ is quite wide, assuming values from 0.62rbto 0.83rb. Introduction of a hysteresis range gives rise, in Figure 311, 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. 311 (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 310, to convex. For a higher meniscus, h/ = 1.5, for the same 4)^ and Ap shown in Figure 313, the radius variation is also not very substantial in magnitude and the waveform is not identical to h<.(t), indicating the nonlinearity that is introduced by <^(t). When hysteresis is present, and for the case of fairly large hysteresis, and tall meniscus. Figure 314 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 314 (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 quasiequilibrium 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 quasiequilibrium situation, and the full fluiddynamic
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 fourthorder 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 nonextrema 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 nonequilibrium 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. Quasiequilibrium 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 31 Schematic of meniscus corresponding to edgedefmed fibre growth process.
PAGE 97
1 Â« . 80 Figure 32 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 32 continued
PAGE 99
82 Figure 33 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 34 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 35 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 36 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 37 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 38 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 312 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 314 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 315 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 oftused quasistationary 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 nondimensional 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 steadystate 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 Onedimensional 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 semiinfmite 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 onedimensional transient phasechange problems.
PAGE 115
98 This method provides a relatively straightforward and simple approach for approximate analysis of onedimensional transient phasechange 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 onedimensional 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 twodimensional geometries to arbitrary shaped, doubly connected twodimensional 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 phasechange 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 fixedgrid approach and the transformedgrid approach. The fixedgrid 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 transformedgrid 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 enthalpyporosity 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 twodimensional convective melting of pure gallium in a rectangular cavity. The enthalpyporosity technique was validated by comparing its results with the experimental results obtained by Gau and Viskanta (1986). A review of available implicit difference enthalpyporosity 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 solidliquid interface and, consequently, appear to have been the first investigators to apply it to a multidimensional convectiondominated 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 quasisteady 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 quasisteady 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 timemarching technique using FEM but was not applied to the moving solidliquid 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 quasistationary 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 singlevalidness 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 slowmoving interfaces. 4.3.1 Assessment of the Quasistationary 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. 41) confined to a
PAGE 122
105 semiinfinite 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 nondimensionalizing 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 secondorder central differences are used for all spatial derivatives. The quasistationary 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 onedimensional 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 secondorder 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) Quasistationary treatment: By invoking the quasistationary 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. 42 (a), both full and quasistationary 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 quasistationary 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 quasistationary 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 = il,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 multiplevalued 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 nondimensional 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. 43, 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. 42. 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 Timestepping 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 CrankNicolson schemes for solving onedimensional 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 41 (a) and 41 (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 41 (a) the accuracy of the CrankNicolson 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 firstorder onesided scheme is used to estimate the thermal gradients at the interface. Consequently, the marker velocity is firstorder 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 CrankNicolson, 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 11 and 12, 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 42, 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 solidliquid interface, the Rayleigh number is 0(10'*) for sapphire, according to the property values listed in Table 42. 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 42 and a surface tension of 0.06 dyne/cmK 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 GibbsThompson effect or the anisotropic surface energy (Pelce 1988).
PAGE 133
116 We are interested in the axisymmetric situation, i.e. rz 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 44. Table 42 gives the physical dimensions and the material properties of sapphire, while Table 43 summarizes the important parameters and their numerical values encountered in the present model based on Table 42. 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 rz 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 (423) 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 kfIl pMu^\ (430 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) <^^ (431) 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 steadystate. 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.[Lr_(/?_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 pseudosource terms, we move them to the right hand side of the equation: where the spatial derivative terms are treated based on the secondorder central difference scheme in the context of secondorder 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 powerlaw 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 nondimensional 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 45. 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 nondimensional operating parameters are given in Tables 42 and 43. 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 YMP 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 46, and the corresponding steadystate solution is shown in Figure 47 and 48,
PAGE 144
127 respectively. Density and thermal diffusivity are taken to be the same for both the melt and solid phases. Figure 47 shows the isotherms of the steadystate 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 steadystate isotherm distribution of St = 0.024 are shown in Figure 48; 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 nondimensionalize 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 44 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 nondimensional speed of interface movement is scaled with St. Accordingly, in nondimensional 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 nondimensional 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 nondimensional 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 49 compares the paths of the computing histories with both scaling procedures from an identical initial condition toward the steadystate 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 nondimensional 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 nondimensional 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 steadystate. As a check, the converged steadystate 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 steadystate 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 quasistationary 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 singlevalue 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 41 Assessment of time stepping schemes for phase change problem. (a) Comparisons of the predicted interface location at r = 0.2 using backward Euler and CrankNicolson 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 CrankNicolson Backward Euler CrankNicolson 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 CrankNicolson 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 CrankNicolson Backward Euler CrankNicolson 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 42 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/cmK) 0.1 0.1 0.1 p, (g/cm') 3.05 3.05 3.05 Cp, (J/gK) 1.26 1.26 1.26 e, 0.9 0.9 0.9 k, (W/cmK) 0.1 0.1 0.1 p, (gW) 4.00 3.05 3.05 Cp, (J/gK) 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/cm2K) 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/cms) 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 43 Definition and magnitude of key dimensioniess parameters based on values given in Table 41 (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 44 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 44 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 41 Exact solution for melting of a solid confined region. to a semiinfinite
PAGE 155
138 Figure 42 Comparison of interface trajectories predicted by the coupled (full equations) and decoupled (quasistationary equations) approaches with three values of the Stefan numbers.
PAGE 156
139 Figure 43 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 44 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 45 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 46 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 47 The grid system and isothermal contours of the steadvstate solution for Case 1: St = 1. '
PAGE 161
144 (a) grid system (b) isotherm (A0 = 0.2) Figure 48 The grid system and isothermal contours of the steadystate 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 meniscusdefined 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 steadystate solution for St = 1 shown in Figure 47 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 sawtoothed 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 51 and 52, 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 threeharmonic 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 51 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 timeaveraged mean values are different from the steadystate values. This phenomenon indicates that there are accumulated effects caused by the nonlinear
PAGE 167
150 physics intrinsic to the process. Finally, Figure 51 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 53 for both the singleharmonic 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 singleharmonic perturbation (Q=500At) exhibits more interface movement than that with threeharmonic 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 multidimensional 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 53. 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 54. Figure 54 (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 54 (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 (28 cm/min) is one to two orders of magnitude larger than that of conventional growth methods, e.g. Czochraiski (Zulehner 1983) or floatingzone (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 steadystate condition represented in Figure 47 is used as the initial condition to start the computations. Two types of perturbation are considered, namely, the sinusoidal and the sawtoothed 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 singleharmonic 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 55. 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 sawtoothed 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 56. Figure 56(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 56(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 54(a) and 56(a). The variations of
PAGE 171
154 H, and with respect to the forcing frequency are presented in Figure 56(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 53(b), Figure 56(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 52. Figure 57 shows the time histories of Up , H, and R, for this case. While Up in Figure 57 and 6^ in Figure 52 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 sawtoothed 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 58. 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 stateoftheart 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 43. 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 steadystate solution. Figure 58, as the initial condition. Figure 59 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 steadystale values. Figure 510 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 nondimensional time scale, which is of order 1 according to choice 2 given in Table 43, 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 510. This phenomenon has also been observed in the actual growth system. Figure 511 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 512, 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 512, 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 513 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 timeaveraged 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 timeaveraged 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 48, 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 513. 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 513. 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 514 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 quasiperiodic 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 515 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 timeaveraged values of and are different from the steadystate 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 51 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 52 Time histories of and R, subject to a threeharmonic 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 singiehannonic permrbaiion in 6^ with Qj = 500 At 0.21 tÂ»l300dt \. t1600dt Â• tÂ» 1700 dt "Â•Â• Â• ... i1800dt n/'S... 0 0.1 Oa OJ 0.4 OJ 0.6 0.7 0.8 (b) a threeharmonic pemirbauon in 6^ with Q, = 500 At, ~ 292 At and Q, = 108 At . Figure 53 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 55 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 56 Characteristics of interface movement with respect to different H ani uTT: ^= "'^''^ P*' between o/ tan
PAGE 186
169 ure 57 Time histories H, and subject to a threeharmonic 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 58 J^^^ories of sawtoot^^^ 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 59 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 513 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 514 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 (nondimensional time steps) (a) Percentage variation of 12 10 8 4 50 100 150 200 250 300 350 400 450 500 Period (nondimensional time steps) (b) Percentage variation of R, Figure 515 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 fluiddynamic consideration is required to assess the implication of such behavior. A new generalized interface tracking technique is developed, which overcomes restrictions of singlevaluedness 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 timeaveraged values of and R, are different from the steadystate 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 BoundaryMoving Mesh Analysis of Phase Change Using Finite Element Method," International Journal for Numerical Methods in Engineering, 23, pp. 591607, 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. N0001490C0060, 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. 75150, Academic Press, New York, 1964. Beckermann, C. and Viskanta, R., "DoubleDiffusive Convection Due to Melting," International Journal of Heat Transfer and Mass Transfer, 31, pp. 20772089, 1988. Bennon, W.D. and Incropera, F. P., "A Continuum Method for Momentum, Heat and Species Transport in Binary SolidLiquid Phase Change SystemsI. Model Formulation," International Journal of Heat and Mass Transfer, 30, pp. 21612170, 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. 150172, Oxford University Press, London, 1975. Brent, A.D., Voller, V. R. and Reid, K. J., "EnthalpyPorosity Technique for Modeling ConvectionDiffusion Phase Change: Application to the Melting of a Pure Metal," Numerical Heat Transfer, 13, pp. 297318, 1988. 183
PAGE 201
184 Brown, R. A., "Theory of Transport Processes in Single Crystal Growth from the Melt," AIChE Journal, 34, pp. 881911, 1988. Carslaw, H.S. and Jaeger, J. C, Conduction of Heat in Solids . London, Oxford University Press, London, 1959. Chemousko, F. L., "Solution of NonLinear Heat Conduction Problems in Media with Phase Change," International Chemical Engineering, 10, pp. 4248, 1970. Clyne, T. W., "Numerical Treatment of Rapid Solidification," Metallurgical Transaction B 15B, pp. 369381, 1984. Correa, S. M. and W. Shyy, "Computational Models and Methods for Continuous Gaseous Turbulent Combustion," Progress in Energy and Combustion Science, 13, pp. 249292. 1987. Crank, J. . Free and Moving Boundary Problems . Oxford University Press, London, 1984. Crank, J. and Gupta, R. S., "Isotherm Migration Method in TwoDimensions," International Journal of Heat and Mass Transfer, 18, pp. 11011106, 1975. Crochet, M. J.,Geyling, F. T. and Schaftingen, J. J. van, "Numerical Simulation of the Horizontal Bridgman Growth. Part 1: TwoDimensional Flow," International Journal for Numerical Methods in Fluids, 7, pp. 2749, 1987. Crowley, A. B., "Mathematical Modeling of Heat Flow in Czochralski Crystal Growth," IMA Journal of Applied Mathematics, 30, pp. 173189, 1983. Cuvelier, C. and Schulkes, R. M. S. M., "Some Numerical Methods for the Computation of Capillary Free Boundaries Governed by the NavierStokes Equations," SIAM Review, 32, pp. 355423, 1990. Dantzig, J. A., "Modeling LiquidSolid Phase Change with Melt Convection," International Journal for Numerical Methods in Engineering, 28, pp. 17691785, 1989. de Gennes, P. G., "Wetting: Statics and Dynamics," Reviews of Modem Physics, 57, pp. 827863. 1985. Derby, J. J. and Brown, R. A., "ThermalCapillary Analysis of Czochralski and Liquid Encapsulated Czochralski Crystal Growth," Journal of Crystal Growth, 74, pp. 605624, 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: ThreeDimensional Flow," International Journal for Numerical Methods in Fluids, 7, pp. 4967, 1987. Duranceau, J. D. and Brown, R.A., "ThermalCapillary Analysis of SmallScale Floating Zones: Steadystate 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. 371400, 1979. Dyson, D. C, "The Energy Principle in the Stability of Interfaces," Progress in Surface and Membrane Science, 12, pp. 479564, 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 FreezingMelting Interfaces in Fluid Flow, "Annual Review of Fluid Mechanics, 15, pp. 293319, 1983. Ettouney H. M. and Brown, R. A., "FiniteElement Methods for Steady Solidification Problems," Journal of Computational Physics, 49, pp. 118150,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. 323340, 1989. Gau, C. and Viskanta, R., "Melting and Solidification of Pure Metal on a Vertical WaU," Journal of Heat Transfer, 108, pp. 174181,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. 647668, 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. 199217, 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. 335342, 1958. Haley, P. J. and Miskis, J. M., "The Effect of the Contact Line on Droplet Spreading," Journal of Fluid Mechanics, 223, pp. 5781, 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. 883891, 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. 1138, 1977. Huh, C.,and Scriven, L. E., "Shapes of Axisymmetric Fluid Interfaces of Unbounded Extent," Journal of Colloid and Interface Science, 30, pp. 323331, 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. 123136, Hemisphere, New York, 1991. Jansons, K. M., "Moving Contact Lines at Nonzero Capillary Numbers," Journal of Fluid Mechanics, 167, pp. 393407, 1986. Joanny, J. F. and de Gennes, P. G., " A Model for Contact Angle Hysteresis," Journal of Chemical Physics, 81, pp. 552562, 1984. Johnson, R. E., "Conflicts between Gibbsian Thermodynamics and Recent Treatments of Interfacial Energies in SolidLiquidVapor Systems," Journal of Physical Chemistry, 63, pp. 16551659, 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. 31632, 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. 289295, 1990. Koplik, J.,Banavar, J. and Willemsen, J. P., " Molecular Dynamics of Fluid Flow at Solid Surfaces," Physics of Fluids, A 1, pp. 781794, 1989. LaBeile, H. E. Jr., "EFG, the Invention and Application to Sapphire Growth," Journal of Crystal Growth, 50, pp. 817, 1980. Lacroix, M., "Computation of Heat Transfer During Melting of a Pure Substance from an Isothermal Wall," Numerical Heat Transfer, B, 15, pp. 191210, 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. 2541, 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. 5998, 1979. Levi, C. G. and Mehrabian, R., "Heat Flow during Rapid Solidification of Undercooled Metal Droplets," Metallurgical Transaction, A 13 A, pp. 221234, 1982. Lightfoot, N. M. H., "The Solidification of Molten Steel," London Mathematical Society Proceedings, 31, pp. 97116, 1929. Lobar, B. L. and Jain, P. C, "Variable Mesh Cubic Spline Technique for NWave 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. 387441, 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. 18511863, 1979. Michel, D. H., "Meniscus Stability," Annual Review of Fluid Mechanics, 13, pp. 189215, 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. 237253, 1971. Neumann, A. V. and Good, R. J., "Thermodynamics of Contact Angles, I. Heterogeneous Solid Surfaces," Journal of Colloid and Interface Science, 38, pp. 341358, 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. 2740, 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. 568581, 1977. Ostrach, S., "Fluid Mechanics in Crystal Growth 1982 Freeman Scholar Lecture," ASME Journal of Fluid Engineering, 105, pp. 520, 1983. Ostrach, S., "LowGravity Fluid Flows, " Annual Review of Fluid Mechanics, 14, pp. 313345, 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. 370377, 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. 641651, 1976. Roache, P. J., "On Artificial Viscosity;" Journal of Computational Physics, 10, pp. 169184, 1972.
PAGE 206
189 Robison, A. C, "A Survey of Optimal Control of DistributedParameter Systems," Automatics, 7, pp. 371388, 1971. Roe, P. L., "CharacteristicBased Schemes for the Euler Equations," Annual Review of Fluid Mechanics, 18, pp. 337365, 1986. Sachs, E. M., "Thermal Sensitivity and Stability of EFG Silicon Ribbon Growth," Journal of Crystal Growth, 50, pp. 102113, 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. 453492, 1989. Saitoh, T., "Numerical Method for MultiDimensional Freezing Problems in Arbitrary Domains," ASME Journal of Heat Transfer, 100, pp. 294299, 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. 445473. Schwartz, A. M. and Tejada, S. B., "Studies of Dynamic Contact Angles on Solids," Journal of Colloid and Interface Science, 38, pp. 359375, 1982. Shingu, P. H. and Ozaki, "Solidification Rate in Rapid Conduction Cooling," Metallurgical Transaction A 6A, pp. 3337, 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. 11571174,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. 241271, 1990. Shyy,W., "A Study of Finite Difference Approximations to SteadyState, ConvectionDominated Flow Problems," Journal of Computational Physics, 57, pp. 415438, 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. 303308, 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. 247261, 1991. Shyy, W. and Chen, M.H., "SteadyState Natural Convection with Phase Change," International Journal of Heat and Mass Transfer, 33, pp. 25452563, 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. 12291245, 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 922902, 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. 287313, 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. 520526, 1977. Sullivan, J. M. Jr. and Lynch, D. R., "NonLinear Simulation of Dendritic Solidification of an Undercooled Melt, " International Journal for Numerical Methods in Engineering, 25, pp. 415444, 1988. Surek, T., Coriell, S. R. and Chalmers, B., "The Growth of Shaped Crystals from the Melt," Journal of Crystal Growth, 50, pp. 2132, 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. 255279, 1975. Tatarchenko, V. A. and Brenner, E. A., "Crystallisation Stability during Capillarv Shaping, I and II, "Journal of Crystal Growth, 50, pp. 3350, 1980. Thomas, P. D.,Ettouney, H. M. and Brown, R. A., "A thermalCapillary Mechanism for a Growth Rate Limit in EdgeDefined FilmFed Growth Rate Limit Silicon Sheets," Journal of Crystal Growth. 76, pp. 339351, 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. 10071019, 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. 117, 1992. Ungar, L.H. and Brown, R. A., "Cellular Interface Morphologies in Directional Solidification. The OneSided Model," Physical Review B, 29, pp. 13671380, 1984. Viskanta, R., "Heat Transfer During Melting and Solidification of Metals," ASME Journal of Heat Transfer, 110, pp. 12051219,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. 845877, 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. 153222, 1983. VoUer, V.R., "Implicit FiniteDifference Solutions of the Enthalpy Formulation of Stefan Problems," IMA Journal of Numerical Analysis, 5, pp. 201214, 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. 545556, 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. 875898, 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. 141153, 1992. Warming, R. and Hyett, B., "The Modified Equation Approach to the Stability and Accuracy Analysis of FiniteDifference Methods," Journal of Computational Physics, 14, 159179, 1974. Warming, R. F.,Kutler, P. and Lomax, H., "Secondand ThirdOrder Noncentral Difference Schemes for Nonlinear Hyperbolic Equations," AIAA Journal, 11, pp. 189196, 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. 131153, 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. 327356, 1987. Zulehner, W., "Czochralski Growth of Silicon," Journal of Crystal Growth, 65, pp. 189213, 1983.
PAGE 210
BIOGRAPHICAL SKETCH ShinJye 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. ChenChi 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 UTF8
REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchemainstance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd
INGEST IEID ER828U6WP_OKN3X8 INGEST_TIME 20150422T21:51:45Z PACKAGE AA00029996_00001
AGREEMENT_INFO ACCOUNT UF PROJECT UFDC
FILES
. Â— ^ CO il Â« o u ^
PAGE 70
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 infmite extent. On the other hand, in the crystalpulling 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
PAGE 71
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 (f>^ 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 4>^ is greater than that formed by receding liquid, say (f)^. Thus, the apparent angle assumed under static condition by the liquidgas surface at the trij unction is either or <^r . 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 to <^R or viceversa, depending on the direction of impending motion.
PAGE 72
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 f^ 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 LaplaceYoung equation for the equilibrium meniscus profiles. The boundary conditions are chosen to correspond to that of the edgedefined filmfed 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.
PAGE 73
56 3.2 Configuration. Assumptions and Formulation The configuration chosen here is shown in Figure 31, where hÂ„ re, and rt 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 LaplaceYoung equation for the liquidgas interface is applied, which in dimensional form is (Surek et al. 1980, Pitts 1976, and Huh and Scriven 1960) = Ap5z Ap where f is the radial coordinate defining the surface of the meniscus, is the surface tension of the liquidgas interface, Ap is the difference in density between liquid and gas phases, and Ap is the difference in applied pressure between the liquid, pi , and gas phases, Pg , i.e.,Ap = Pi 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
PAGE 74
57 a= JJl. (3.2) where, for example, for sapphire (AI2O3), 7,g = 0.69 N/m, Ap = 3.4x10? kg/m^ , and hence the capillary length is about 4x10" m. the bond number, which relates radius of base or die to capillary length, and controls relative importance gravity forces is given by , bo="^" (33) a' we specify aspect ratio, hjty, bo, each case. meniscus height calculated from ratio r^. nondimensionalize all physical quantities as follows,="^'s" representing nondimensional values, z*="z/a" ,r="f/a" ap*="Ap/(Apga)," ap '="1." again, using sapphire an example, reference pressure p,="Ap" g a about 130 n m^ thus, nondimensionalized profile equation, after having dropped *'s is, fzz 1 t r="^^p" (3.4) this equation seemingly invariant with respect since number does not explicitly appear eq. (3.4). misleading impression, influence profiles via term z, term. for h, rb, effect

