NASA Contractor Report 3962

An Exact Solution for the

Solidification of a Liquid
Slab of Binary Mixture

B. N. Antar, F. G. Collins,
and A. E. Aumalis
The University of Tennessee Space Institute
Tulahoma, Tennessee

Prepared for
George C. Marshall Space Flight Center
under Contract NAS8-34268

National Aeronautics
and Space Administration
Scientific and Technical
Information Branch

1986
 FOREWORD
The work reported here concerning formulation and solution
of a two phase mixed Stephan problem is relevant to alloy
solidification and crystal growth processes such as those being
investigated in low gravity experiments aboard orbiting
laboratories. The work was supported by the Atmospheric Science
Division of the Systems Dynamics Laboratory at the Marshall Space
Flight Center.
 TABLE OF CONTENTS

CHAPTER

I . .......................
INTRODUCTION 1

I1 . PROBLEM FORMULATION . . . . . . . . . . . . . . . . . . . 4

The Solutions of Neumann and Stefan . . . . . . . . . . 9

Numerical Technique. the Method of Meyer . . . . . . . . 12

Meyer's Method for the Mixed Problem . . . . . . . . . . 17

Calculation of Initial Condition . . . . . . . . . . . . 19

I11 . SOLUTION AND RESULTS . . . . . . . . . . . . . . . . . . . 25

Numerical Method . . . . . . . . . . . . . . . . . . . . 25

Analytical Method . . . . . . . . . . . . . . . . . . . 28

Computation Results . . . . . . . . . . . . . . . . . . 37

LIST OF REFERENCES . . . . . . . . . . . . . . . . . . . . . . . 68

. . . .

iv
 LIST OF FIGURES

FIGURE PAGE

1. Sketch of theoretical model................ 5

2. Pb-Sn Phase diagram .................... 8

3. Comparison of interface position vs. time for numerical

and analytical techniques. Solid = analytical,

Circles = numerical. .................. 39

4. Comparison of temperature vs. time for numerical and

analytical techniques. Solid = analytical, Circles =

numerical. ....................... 40

5. Comparison of concentration vs. time for numerical and

analytical techniques. Solid = analytical, Circles =

numerical. ....................... 41

6. Concentration distributions at P = 10
S
-4
......... 42

7. Concentration distributions at P = 10
S
-5
......... 43

8. Concentration distributions at P = 10
S
-6
......... 44

9. Interface vs. time at P =

S
.............. 45

10. Concentration vs. time at P = 10

S
-4
............ 46

11. Temperature vs. time at PS = 10

-4
............. 47

12. Temperature vs. interface at PS = 10

-4
.......... 48

13. Interface vs. time at PS = .............. 49

14. Concentration vs. time at P = 10

S
-5
............ 50

15. Temperature vs. time at PS = 10

-5
. . . . . .-... . . . . . . 51

16. Temperature vs. interface at PS = 10

-5
........... 52

17. Interface vs. time at P = 10

S
-6
.............. 53

V
FIGURE PAGE

18. Concentration vs. time at P

S
= 10
-6
............ 54

19. Temperature vs. time at P

S
= 10
-6
............. 55

20. Temperature vs. interface at P = 10-6 . . . . . . . . . . 56

S

21. Comparison of interface vs. time for R = 1.00 and

R=1.06 ........................ 58

22. Comparison of temperature vs. time for R = 1.00 and

R=1.06 ........................ 59

23. Comparison of concentration vs. time for R = 1.00 and

R=1.06 ........................ 60

24. Temperature and concentration distributions for time =

25 - 100 with Ps = ................. 61

25. Temperature and concentration distributions for time =

125 - 200 with Ps = ................ 62

26. Temperature and concentration distributions for time =

25 - 100 with Ps = ................. 63

27. Temperature and concentration distributions for time =

125 - 200 with Ps = 10

-5
................ 64

28. Temperature and concentration distributions for time =

25 - 100 with Ps = ................. 65

29. Temperature and concentration distributions for time =

125 - 200 with Ps = ................ 66

-6
30. Comparison of analytical method at P S = 10 with
intuitive results from reference [ 1 3 ] .......... 67

vi
 LIST OF SYMBOLS

depth of slab

d solute diffusivity

C concentration

C specific heat
P
h heat transfer coefficient

k thermal conductivity

i; partition coefficient

11 latent heat

slope of solidus-liquidus line in phase diagram

interface

temperature

t time

z distance

IC thermal diffusivity

P density

vii
 CHAPTER I

INTRODUCTION

Heat transfer problems dealing with the melting or freezing of

a substance have received a lot of attention for more than a century.

The problem is characterized by a moving surface of separation between

the two phases, with release or absorption of latent heat at this inter-

face. Of interest are the position of this interface with respect to

time and the temperature distributions in both phases. This class of

problems are generally referred to as the Stefan problem.

The simplest problem in this area is the one-dimensional solid-

ification in a semi-infinite region. The first published work is by

G. Lame and B. P. Clapeyron (1831) [17],l which dealt with determining

the thickness of the solid crust generated by the cooling of a liquid

initially at constant temperature. They found the position of the

interface to move proportionally to the square root of time, but they

did not find the constant of proportionality. The first known solu-

tions for this problem were found by Neumann (1860) and Stefan (1889)

[ 6 ] . The solution takes the form of error integrals and the position

of the interface was found to be S(t) = 2 X G , where K~ is the thermal

diffusivity, and X is determined by the conservation of energy equation

at the interface.

Further complications arise when the material freezing is an

alloy. The classical Stefan problem normally has three unknowns,

'Numbers in brackets refer to similarly numbered references in

the List of References.

1
namely, the interface position and the temperatures of both phases.

The melting temperature of an alloy is strongly dependant on the com-

position at the interface. Therefore, freezing (or melting) can occur

over a range of temperature. That, and the concentration distributions

of the two phases, adds three more unknowns to the problem. Work done

by Tao [ 1 8 ] considers this problem with an arbitrarily prescribed heat

flux at the lower boundary. His closed form solution is comprised of

an infinite sum of error functions.

When considering a slab of finite depth, the similarity trans-

formation, in general, cannot be applied and no known closed form

solution exists. Boley [ 3 ] has formulated a solution for the system

as it first undergoes solidification. However, once the effects of

freezing have reached the upper surface the solution is no longer

valid. To describe the system for all time, we must resort to numer-

ical approximation techniques. Meyer [11,12] has developed an

algorithm using the Implicit Imbedding technique that solves the prob-

lem of the slab, but does not take into consideration the case of alloy

solidification.

The study presented here considers alloy solidification in a

finite slab, with heat dissipation from both the upper and lower sur-

faces by convection. Though the governing equations are linear, these

interfacial conditions are not, making this a non-linear problem.

Another complication will be added, that of allowing for density dif-

ferences between solid and liquid phases. This causes a moving upper

boundary as well as the moving interface. We see this solution as an

extension of Meyer's method. Two methods of solutions will be

2
presented. One, a purely numerical technique, makes use of a Runge-

Kutta integrator to solve the equations generated by Meyer's method.

This solution failed at the low diffusion rates of solids because of

the stiffness of the equations. The second solution is an analytical

method. Closed form solutions were found of the basic equations,

though quadrature techniques were still required. This enabled us to

reach lower diffusion rates, though not as low as desired, because of

accuracy limitations of the quadrature technique.

Stefan type problems have significant importance to many indus-

trial processes. In the production of steel billets [ 1 4 ] , carbon

within the metal oxidizes at the surface setting up strong concentra-

tion gradients that cause more carbon to diffuse to the surface. This

outer layer is no longer of the same type steel and would have to be

ground off. Determining and minimizing the depth of this layer is

necessary for the reduction of costs. Corrosion problems are similar

in behavior. In the manufacture of glass [ 1 4 ] , heat transfer rates

must be optimized to limit crystal growth and reduce the formation of

gas bubbles. Recently, much interest has been generated in the

possibility of processing materials in space. In the reduced gravity

environment, it would be possible to produce large crystals of uniform

properties and to manufacture materials with unique properties.

3
 CHAPTER I1

PROBLEM FORMULATION

I n t h i s c h a p t e r w e w i l l d i s c u s s t h e f o r m u l a t i o n of t h e b a s i c

e q u a t i o n s f o r t h e two phase mixed S t e f a n problem. Exact s o l u t i o n s by

Neumann and S t e f a n and t h e numerical t e c h n i q u e by Meyer w i l l be

reviewed. For t h e mixed S t e f a n problem w e w i l l c o n s i d e r a s l a b con-

s i s t i n g of a b i n a r y l i q u i d a l l o y , i n i t i a l l y a t c o n s t a n t t e m p e r a t u r e ,

The s l a b is of f i n i t e d e p t h , do, and i n f i n i t e width ( F i g u r e 1).

To.
Heat i s t o b e removed from b o t h t h e upper and lower s u r f a c e s , w i t h

s o l i d i f i c a t i o n s t a r t i n g a t t h e lower s u r f a c e . The problem i s charac-

t e r i z e d by a moving i n t e r f a c e , S ( t ) , between t h e s o l i d and l i q u i d

phases, w i t h l a t e n t h e a t release a t t h e i n t e r f a c e . The e q u a t i o n s w i l l

a l l o w f o r d i f f e r e n t d e n s i t i e s f o r t h e s o l i d and l i q u i d c a u s i n g a

volume change and t h e r e f o r e a moving upper s u r f a c e d u r i n g f r e e z i n g .

The upper and lower s u r f a c e s w i l l b e c l o s e d t o any mass flow and mass

d i f f u s i o n i n t h e s o l i d and l i q u i d w i l l b e i n v e s t i g a t e d . S i n c e t h e sub-

s t a n c e i s an a l l o y , t h e f r e e z i n g temperature w i l l b e dependant on t h e

composition a t t h e i n t e r f a c e . A s k e t c h of t h i s model is shown i n

F i g u r e 1.

The governing e q u a t i o n s f o r t h e temperature and c o n c e n t r a t i o n

d i s t r i b u t i o n s when t h e r e i s a change i n volume d u r i n g s o l i d i f i c a t i o n

(ps/p # 1) t a k e t h e form [61:

R

4
 0
QI
5

0 n
a u
W II
II VI
N
N

5
in the liquid, and

2
--
aTS a Ts
-
at
- K
s 2 (3)
az

2
-- a cs
at - ds 2 (4)
az

in the solid. C denotes the solute concentration and T the temperature.

S(t) is the interface position and p , K and d denote the density, thermal

diffusivity and solute diffusivity. The subscripts s and R denote solid

and liquid, respectively.

The interface conditions are:

aTS aTR dS
ks az - kt az = psR dt (5)

acS
Psds Ti-- PRdR az -
Equations (5) and ( 6 ) are the conservation of energy and mass equations
c,

at the interface. k, k and R are the thermal conductivity, the partition

coefficient and the latent heat. Also at the interface, both the solid

and liquid must be at the equilibrium melting temperature, Tm, and the

proper value of solute concentration which may be described by

6
 Tm = Ti - mC (7)
II

Equation (7) implies that the melting temperature of the alloy

is a linear function of concentration. Looking at the phase diagram,

Figure 2 [ l ] , we can see that this is a simplification of the actual

relationship. The decrease in the melting temperature with increase

in concentration can be rather closely approximated by a linear func-

tion. However, at a certain concentration the melting temperature

will begin to increase. Since the method of solution does not depend

on this equation, a more complicated expression could be substituted,

but will not be for this work.

The boundary conditions at the upper and lower surfaces take

the form:

atz=O

aT = h2 (T,
-
- k g az
- T2)

7
 t
v)

0
cn

0
a
L
i=I .-c
I-
c

U Q)o
0 FF
c
0
t
Q) C c
I Q)
? os El
rn
I
$0 O & P
(D n b

E .-0
E N
v) 0
c a,
I 40
n tn
n
0
d-

%
0
N

0
a
n

8
where subscripts 1 and 2 denote lower and upper surfaces, respectively.

Equations (9) and (11) represent heat dissipation from the surfaces by

convection into a medium at temperatures T and T2, respectively.

1

The Solutions of Neumann and Stefan

The first solutions for the Stefan problem, in a semi-infinite

space, were found by Neumann in 1860 and Stefan in 1889 [ 6 ] . The

governing equations are (l), (2), and interface condition (5). They do

not allow for a density change between the solid and liquid, i.e.,

Ps
-- Pg’ The boundary conditions at zero and are

Ts = 0 at z = 0

The solutions

TR = T - Berfc (- = >
2 f i
S

Ts = Aerf (- = >
2 f i
S

satisfy the basic equations and boundary conditions at zero and my

where A and B are constants of integration. Using the interface con-

dition

9
 Ts=T = T a t z = s
R m

we have

S S
Aerf (- ) = T - Berfc (-) = Tm
2Q 2 4 7

For t h i s t o b e v a l i d f o r a l l t , w e must have

S
-= X = constant
2 q -

or

s = 2hQ

r e w r i t i n g t h e s o l u t i o n a t t h e i n t e r f a c e , w e have

Aerf(X) = T - Berfc ( A
p)
KR = Tm

which d e t e r m i n e s t h e c o n s t a n t s A and B as follows:

m
A =
erf (A)

T - Tm
B -

erfc
R

Now, u s i n g t h e i n t e r f a c e c o n d i t i o n ( 5 ) , w e can e v a l u a t e t h e c o n s t a n t A :

10
 2
exp(- -
K
R
exp(-A ) - -k T
m =- ARJ;;
C T
e r f (A) m
S ps

S t e f a n ' s s o l u t i o n assumes t h e l i q u i d i s i n i t i a l l y a t t h e m e l t i n g

temperature (T = Tm), g i v i n g t h e e q u a t i o n f o r A a s

2
Aexp(A ) e r f ( A ) = - u

a?G

Adding t h e change i n d e n s i t y t e r m (p / p
s R
- l), so that the basic

e q u a t i o n matches e q u a t i o n , ( l ) , changes t h e s o l u t i o n i n t h e l i q u i d

r e g i o n t o t h e following:

Ts = T - B e r f c [-
2J"R"
+ A ( - -p~
PR
1) J21 R

with

T - T
m
B =

er f c (A-

and t h e e q u a t i o n t o determine A :

11
 Numerical Technique, the Method of Meyer

When considering a slab of finite depth, a numerical technique

must be employed. In a series of papers Meyer [11,12] has developed

one such technique using the Invariant Imbedding Method. This tech-

nique also allows for the use of non-linear boundary conditions. The

governing equations are (1) and ( 3 ) , with interface conditions (5),

and boundary conditions (9) and (11). Since numerical methods are

used, it is convenient to first nondimensionalize the equations. The

following substitutions were made:

c = -C Lewis Number
cO

cO
M TI - To
m = -ds
= Lewis Number
ps K
S

K
K = KS hido
B. = -= Biot Modulus
R = kR

12
 z = -z
do

k T2 - To
S
A=- Y =
kg T1 - To

T - To K
o =
T1 - To

After nondimensionalization, the governing equations, interface

and boundary conditions take the form

aos a 2os
-- - K-
aT 2
aZ

aos
A---=
aoR RL -
dS
az az dT

Meyer's method consists of first discretizing the time operator

in the following manner:

13
yielding a set of ordinary differential equations, solvable at discrete

time steps, N. Next, equations (13) and (14) are rewritten as a first

order system.

N
0% - oIIN- 1 - - N- 1
- (R 1) (SN S ) uII = u’n.
AT AT

osN - osN- 1
AT = KU’, (19)

where prime denotes differentiation with respect to z . With the con-

servation of energy equation at the interface taking the form

Now, the Invariant Imbedding technique is used to solve for U

R
and Us. The solution of U is assumed to take the form
R

N
+ zII

then

14
 N'

N N N'
*'; = Y
;
1 + Y
R + zR

These e q u a t i o n s are s u b s t i t u t e d i n t o (18) t o y i e l d

N- 1
N' N N N
AT AT FRYR OR - FRZQ

where

- (R - 1) (SN - SN-')
FR - AT

rearranging,

N- 1
0
+ FRYRN ORN - - FRZR --
Y '0 N + (YRN ) 2 0; = Z N' YRN Z RN - N R
R R R AT

from t h i s we g e t t h e e q u a t i o n s f o r Y and Z
R R'

Y
N'
= - N 2
(YR ) N
+ -AT1
R -

N-1
ZR
N'
= - N
YE
N N --uR
- FRZR AT

The boundary c o n d i t i o n s are d e r i v e d s i m i l a r l y .

N
UaN
= YII" 0
R
+ Z RN

i s s u b s t i t u t e d i n t o (17) t o y i e l d

15
and theref ore

N
ZR (Zf) = - B2Y

where Zf '= 1 - (R - 1)s

N
. Similarly,

usN =
N
Ys osN + zsN

yields the following equations and boundary conditions:

dY ,"
- - - -- 1 N 2
dz KAT us1

N = -B 1
at z =
0Ys A

N- 1
-dZ=
," - N N OS
dz
Ys zs - ___
KA 'c

After solving equations (21)-(24), the temperature distributions can

be obtained by solving the following equations:

16
 N
-doR
- -
dz Y; 0; + z RN

N
dos
- dz
- - Ys
N
osN + zsN

Meyer's Method for the Mixed Problem

Meyer's method can easily be applied to alloy solidification.

Equations ( 2 ) and ( 4 ) , interface conditions ( 6 ) , ( 7 ) and (8), and

boundary conditions (10) and (12) are nondimensionalized, time operator

discretized, and then written as a first order system yielding:

c; - N- 1
CR
- 1) - N- 1
- (R (SN S
VR = P V '
AT AT R R

Vs = C I S

csN - c;-1
= KP V'
AT s s

Om = Oi - MCR

cs = iCR

17
 ac,"
-= 0 (32)
az

ac,N
-- - 0 (33)
aZ

V and Vs are assumed to have solutions of the form

R

N N N
= XI1 CI1 + RR
N N
V," = Xs Cs + R,"

Then the equations to solve for X R , Xs , RE and R," are

at z = Zf x; = 0

dR:
-- - - [ FR N ] RN ~ --
P + X ~
PRAT (35)
dz
R

at z = Zf :R = 0

at z = 0 x,"=o

18
 d
R
: cSN-1
--
dz
-- xsN N
Rs
--U S A ~ (37)

at z = 0 R:=O

After solving equations ( 3 4 ) - ( 3 7 ) , the concentration distributions can

be obtained by solving the following equations:

d
C
:
-- - XgN CRN + RgN
dz

d
C
:
-- - XsN CsN + RsN (39)
dz

Calculation of Initial Condition

In order to solve equations ( 2 2 ) , ( 2 4 ) , (35) and (37), temper-

ature and concentration distributions from the previous time step are

needed. T h i s necessitates initial distributions for both the temper-

ature and concentration at the onset of solidification. For the present

case, we will assume that no mass diffusion occurs prior to freezing.

Thus, the concentration distribution is taken to be

throughout the liquid, up until solidification starts.

Since the initial temperature, Oo, may not initially be equal to

the melting temperature, such an assumption may not be made. We

Qm’
must therefore derive an equation for the temperature as a function of

19
p o s i t i o n and t i m e , which i s a p p r o p r i a t e t o t h e s p e c i f i c c o o l i n g a t t h e

upper and lower s u r f a c e s . The governing e q u a t i o n [ 6 ] f o r t h i s i s

w i t h boundary c o n d i t i o n s

aog
-az- - Bl(Oll - 1) at z = 0 (41)

aoR
- = - B2(Og+Y) at z = 1
az

and w i t h i n i t i a l c o n d i t i o n s 0 = 0 a t T = 0. The s o l u t i o n i s d e r i v e d
R
u s i n g s e p a r a t i o n of v a r i a b l e s . The s o l u t i o n of t h e homogenous p a r t

t a k e s t h e form

then

ZT' = Z"T

or

T' Z" 2
-A
T Z

20
 _
dT -
- -A 2
T -dLZ
- - -A
2
z
dT 2
dz

2
T = Alexp(-A T) Z = A2sin(Xz) + A 3c o s ( X z )

and

OR = [A2sin(Az) + A3cos(Az)] exp(-A 2 T )

(43)

aoR
-- - A[A2cos(Az) - A3sin(Az)] exp(-X 2 T ) (44)
32

a t z = O

2 2
AA2exP(-A T) = B 1A 3exp(-A T)

B1
A2 = A A3 (45)

a t z = 1

X[A2cos(A) - A3sin(A)] exp(-A 2 T) = - B [A s i n ( A ) + A3cos(A)] exp(-A 2 T)

2 2

using (45)

21
 lA3
B1A3cos(X) - XA3sin(X) = - B2 [ -
X
sin(A) + A3cos(X) 3

2
X - B1B2
(B1 + B2)cos(X) = [ h sin(X) 3

So the solution of the homogenous part is

Where the eigenvalues,

, are determined from (46).
For the non-homogenous part, assume the solution takes the form

Then, using (41) and (42) we get

B = B1(a - 1)

and

combining and solving for a,

22
 B1(l + B2) - B2Y
a =
B1(l + B2) + B2

Now for t h e t o t a l s o l u t i o n we have

Now u s i n g t h e i n i t i a l c o n d i t i o n

we can s o l v e for t h e
%

0 = a + Bz + 1 p” [x
B1
N
sin(A
N
z) + cos(A,z)]
N= 1

- a - Bz = C %EN(z)
N- 1

[ B 12 + A i + % (A N -B; ) s i n (2AN) + 2B1sin 2 (A,)]

AN

23
NOW, to calculate the distribution for the first state, N = 0, we must

find T = such that the temperature at z = 0, is equal to the melt-

ing temperature,
Om
. The initial temperature distribution is then

equal to O ( Z , T
sol)

24
 CHAPTER I11

SOLUTION AND RESULTS

Numerical Method

In order to solve for the temperature and concentration distri-

butions in both the liquid and solid phases, equations ( 2 1 ) - ( 2 4 ) and

(34)-(37) must first be solved. These equations comprise of a set of

coupled, nonlinear ordinary differential equations with their proper

boundary conditions. Such a set of boundary value problems are most

suitably solved with a Runge-Kutta-Fehlberg [ l o ] integrator. For the

present problem, a RK7(8) was used. Before integrating these equa-

N
tions, a new interface position, S , and concentration,,:C must be

assumed. Once these equations are solved, the temperature and con-

centration gradients at the interface may be obtained from equations

( 2 5 ) , ( 2 6 ) , ( 3 8 ) and ( 3 9 ) . The conservation of energy ( 2 0 ) and mass

( 2 9 ) equations at the interface are then used to check these results.

N
If the initial guess on SN and C was in error, these equations will
R
not be satisfied. A better estimate for the correct values of SN and

:C can be obtained from these equations through a Newton-Rhapson

method. This process of evaluating the temperature and concentration

gradients and reestimating the interface position and concentration

is repeated until convergence to the correct values occurs. Now,

equations ( 2 5 ) , ( 2 6 ) , ( 3 8 ) and ( 3 9 ) can be integrated for the tem-

perature and concentration profiles from the interface to the lower

surface for the solid and from the interface to the upper surface for

25
the liquid. The time is then advanced by one time step, AT, and the

process repeated for the new interface position.

This method (Antar [2]) gave satisfactory results down to a

Lewis number, Ps, of order 10

-4
. At lower Lewis numbers overflow

problems occurred due to the stiffness of equation (36). This can

easily be seen by using the following values as an example:

-5
P = 10
S

K = 2.7

AT = 1

To demonstrate the problem, we will use a fourth order Runge-Kutta

technique as outlined in ref [ 4 ] , and the above values with the

boundary condition of zero,

if

y’ = f(z,y) with a -
< z -
< b

y(a> = boundary value

h = -b - a where N = number of intervals

N

then an approximate solution to y is w, found in the following manner:

w = boundary value
0

kl = hf(xi, wi)
 1
k 2 = hf(xi + -h2 , wi +-k )
2 1

1
k3 = hf(x
i
+ h? , w i +yk2)

k4 = hf(xi+l , wi + k3)

= -
1 (kl + 2k2 + 2k3 + k4)
W i+l 6

f o r each i = 0, . . . ,N - 1.

Applying this to equation (36) with 50 intervals we arrive at the

following values:

For i = 0

2
kl = +7.407 x 10

k2 = - 2 . 0 0 3 x 10 3

kg = - 1 . 9 3 1 x 104

k4 = - 7 . 4 6 0 x 10
6

and

6
w1 = - 1.250 x 10

For i = 1

10
kl = -3.127 x 10

27
 18
k2 = -4.889 x 10

k3 = - 1.195 x 1035

k, = -2.865 x 1068
4

The calculation for k causes overflow on the VAX 11/780, which has an
4
upper bound of order 10
38
.
Analytical Method

Inspection of equations (21) and ( 3 4 ) reveals t h e m to be of

the following general type:

There exists a closed form solution [ 1 6 ] for this type of equation.

Finding closed form solutions avoids the use of numerical integrators,

therefore rendering the problem of stiffness more manageable. The

solution is found by using the transformation

U'
Y=,

The equations then take the form

U" + Bu' - AU = 0

whose solution is found in the following manner.

The characteristic equation is

28
 r2 + Br - A = 0

= -
1 (-B 2 B + 4A)
7
r 1,2 2

u = Clexp(rlz) + C2exp(r2z)

then

rlexp(rlz) + Cr2exp(r2z)
= exp(rlz) + Cexp(r2z)

where

Similarly, equations (23) and ( 3 6 ) have the form

and can be solved using the same transformation

U'
Y =y

yielding

u" = Au

u = clexp(Jji;Z) + C2exp(-fiz)

29
with

Equations ( 2 2 ) , (24-26), ( 3 5 ) , (37-39) are of a more complicated t y p e

-dy- - - F ( z ) y - G(z)
dz

b u t can b e e a s i l y solved u s i n g t h e V a r i a t i o n of P a r a m e t e r s t e c h n i q u e

[91

z Z 2'
y = exp(- F(z)dz) [ C - I G(z')exp (IF ( z ) d z ) d z ' ]
Z Z Z
i i i

Equations (22), ( 2 4 ) , (35) and (37) c o n t a i n f u n c t i o n s whose

v a l u e s are known o n l y a t c e r t a i n p o i n t s , r e q u i r i n g q u a d r a t u r e evalu-

a t i o n s t o determine t h e i r s o l u t i o n s . These t e c h n i q u e s , though, are

known t o b e h i g h l y s t a b l e , a l l o w i n g s o l u t i o n s t o be o b t a i n e d a t lower

L e w i s numbers. The same i t e r a t i v e t e c h n i q u e as d e s c r i b e d earlier i s

used t o determine t h e i n t e r f a c e p o s i t i o n and t h e temperature and con-

c e n t r a t i o n p r o f i l e s a t each t i m e s t e p .

Now t h e s o l u t i o n s of e q u a t i o n s (21)-(24) and (34)-(37) are:

rlexp(rlz) + C 1r 2e x p ( r 2 z )
Y; =
exp(rlz) + Clexp(r2z)

30
 L L

N- 1
= exp(-
Zf
(F,+Y;)dz) [ C2- 1
Zf OR
-
AT
exp( jzf (FR+Y; )dz)dz'] (49)
Z Z 2'

C2 = - B2Y

6K - B,I
c5 =
bK + B1

N- 1
2 0S 2'
N ' N
Zs = e x p ( - I Ys dz) [ C 6 - l - e x p ( l Y;dz)dz' ]
0 0 KAT 0

C6 =
- B1
-A

31
 1 FR
r 1Y2 =T[-P
R
f

N
RR = exp(-
FR
I"f (-+Xx, N
)dz) f
z cN-l
*exp(-
FR
I "f (-+xR N
)dz)dz' (53)
z pR z z' pR

N =- 6z )
tanh( - (54)
S
?T
S

Now t h e s o l u t i o n t o e q u a t i o n s (25), (26), (38) and (39) can b e found

u s i n g t h e V a r i a t i o n of Parameters technique [ 9 ] .

N z
OR ( z ) = e x p ( I' N
YR dz) [ C g + I' N
Z R exp( -f N
YR dz)dz' 3 (56)
z z 2.
i i 1

N ' N N ' N
Os (2) = exp( Ys dz) [ Cl0 +Iz Z s exp( - Ys dz)dz' ] (57)
z. z z
1 i i

z z'
N N
cs (2) =exP( I Xs [ CI2 + I' Rs
N
exp( -I N
Xs d z ) d z ' ] (59)
z z z
i i i

32
 The c o n s t a n t s , C and Cll, r e p r e s e n t t h e temperature and con-
9
c e n t r a t i o n of t h e l i q u i d a t t h e i n t e r f a c e . The c o n s t a n t s , Cl0 and

C12, r e p r e s e n t t h e temperature and c o n c e n t r a t i o n of t h e s o l i d a t t h e

interface. Cg and C must be e q u a l . S i n c e only t h e c o n d i t i o n s a t

10
t h e i n t e r f a c e are known, t h e i n t e g r a t i o n s i n e q u a t i o n (56) and (58)

are performed from t h e i n t e r f a c e up, and, i n e q u a t i o n s (57) and ( 5 9 ) ,

from t h e i n t e r f a c e down.

Asymptotic approximations w e r e used whenever p o s s i b l e t o

reduce t h e chance of overflow e r r o r s . For example, e q u a t i o n (50) i s

e v a l u a t e d as

N
Ys = 6

whenever

Also, exponents were combined a s much a s p o s s i b l e . I n equation (59),

f o r example,

w a s e v a l u a t e d f i r s t , t h e n t r e a t e d as a c o n s t a n t s o t h a t i t could b e

c a r r i e d w i t h i n t h e i n t e g r a l and added t o t h e exponent b e f o r e t h e

exponents are e v a l u a t e d .

Many of t h e i n t e g r a l s could n o t b e solved e x a c t l y because t h e y

c o n t a i n t h e c o n c e n t r a t i o n o r temperature f u n c t i o n s whose v a l u e s are

known o n l y a t c e r t a i n p o i n t s . They were t h e r e f o r e solved u s i n g a 15

33
point Gauss-Legendre quadrature, with maximum interval size of 0.1,

and a cubic spline interpolation [ 4 ] between the known points.

b 15 ri(b - a) +b +a
I F(x)dx = (b - a) C wiF (
2 )
a i=1

The roots and weights (r.,w ) were obtained from Table 2 . 2 in ref
1 i
[51

This method gave good results down to a Lewis number, Ps, of

order Below this value, problems occur due to the accuracy of

the integrator, when applied to equation (59). The error occurs when

is evaluated. To demonstrate this, using the funct-ms

N =- 6 6z
tanh( -)
S
KS

and

tanh(-) 6z

34
we evaluate equation ( 6 0 ) in closed form. Note that these equations

are the same as (54) and (55), except for the missing concentration

term in (55). The exact solution of (60) then is

6
sech(-) - 1
JI;-
S

when evaluated from zero to one. For large values of 6/< the solu-

tion approaches the value of -1.0. Evaluating equation (60) numer-

ically, using a 15 point Gauss-Legendre quadrature with 10 intervals,

and the following parameters

K = 2.7

pS =

AT = 1

z
i
=o

z = l

we get the value of -1.36, which is a difference of .36 from the true

value. Table 1 lists the results of different interval sizing at

different values of P S .

35
 Table 1

The Effects of Different Interval Sizing on the Accuracy

of a Gauss-Legendre Quadrature Technique

Number of Difference
Value Intervals in from
of P, Quadrature Actual Value
-

10 2.33 x

10 3.62 x 10-1

100 1.22

250 1.30

100 3.15 x lo-’

250 5.24

500 4.60 x

Reasonable accuracy is obtained for Ps of order using 250

intervals. However, combined with a 15 point quadrature formula,

this means 3750 function evaluations (actually, equation (55) has an

additional integration, requiring 56,250 evaluations), which means a

high probability of round-off error occurring during the computations.

 Computation Results

As previously discussed, the numerical technique gave good

results down to a Lewis number, Ps, of 10-4 . Alloys having this high

rate of diffusion in its solid phase are rare. For more common
alloys, such as lead-tin (Pb-Sn), the Lewis number is of order 10-10

[7]. Because of this restriction the analytical technique was

derived. With this technique, we were able to get results down to

a Lewis number, Ps, of order 10-6 . The following parameters, based

on the temperature differences

To - T = '5
m

Tm - T1 = '2

T2 - Tm = 4'

were used in the calculations.

K = 2.326 B1 = 1.00

M = - 0.22 B2 = 0.01

L e - 27.7 R = 1.00 or 1.06

-2
PR = 10

k = 0.580 y = - 0.14286

A = 2.703 AT = 1

37
 -4
A more r e a l i s t i c l i q u i d L e w i s number, P k , i s of o r d e r 10 [8]. Compu-

t a t i o n s w e r e performed a t a L e w i s number, of i n o r d e r t o keep

pk ,
a t least t h e r a t i o of l i q u i d t o s o l i d d i f f u s i o n rates r e a l i s t i c . Com-

p a r i s o n of t h e two t e c h n i q u e s a t P = i s shown i n F i g u r e s 3, 4
S

and 5. The i n t e r f a c e p o s i t i o n and t h e temperature and c o n c e n t r a t i o n

a t t h e i n t e r f a c e a s a f u n c t i o n of t i m e are shown t o match e x c e l l e n t l y .

The s o l i d l i n e s r e p r e s e n t t h e a n a l y t i c a l s o l u t i o n w h i l e t h e c i r c l e s

r e p r e s e n t t h e numerical s o l u t i o n .

The n e x t series of p l o t s , F i g u r e s 6 , 7 and 8 , show concentra-

t i o n d i s t r i b u t i o n s a t d i f f e r e n t L e w i s numbers a t t i m e s t e p i n t e r v a l s

of 25, w i t h t h e f i r s t t i m e s t e p a t 2 5 and t h e l a s t a t 200. The i n t e r -

f a c e p o s i t i o n i s seen a s t h e v e r t i c a l l i n e s w i t h t h e s o l i d phase t o

t h e l e f t and t h e l i q u i d phase t o t h e r i g h t . The e f f e c t s of lower

-4
L e w i s numbers a r e v e r y a p p a r e n t . I n F i g u r e 6 , a t L e w i s number 10 ,
t h e r e i s s i g n i f i c a n t d i f f u s i o n i n t h e s o l i d , w h i l e i n F i g u r e 8, a t
-6
L e w i s number 10 , almost no d i f f u s i o n i s n o t i c e a b l e . The o n l y

movement i s a t t h e lower s u r f a c e , which h a s had more t i m e f o r d i f f u -

s i o n t o occur.

F i g u r e s 9-20 show how t h e i n t e r f a c e , c o n c e n t r a t i o n and t e m -

p e r a t u r e v a r y w i t h t i m e and t h e temperature w i t h i n t e r f a c e p o s i t i o n

a t t h e t h r e e L e w i s numbers, and We can see i n

F i g u r e s 9, 13 and 17, t h a t t h e v e l o c i t y of t h e i n t e r f a c e i s decreas-

i n g as t h e a l l o y s o l i d i f i e s . The upward motion of t h e i n t e r f a c e h a s

a l l b u t stopped. T h i s i s due t o e q u a t i o n (7), which d e t e r m i n e s t h e

m e l t i n g temperature of t h e a l l o y . A s seen i n t h e phase diagram, t h e

m e l t i n g temperature d e c r e a s e s , as c o n c e n t r a t i o n i n c r e a s e s , t o a

38
 4

rr
$ 27
aJG
a
rl I1

I-

d I: d
19 c 19

39
 19 a
-Ea
N aC

61
-Ln
4

6
l
-61
4

19
-Ln

-GI
7ll-I- T l - l l -
h
I c rl
a
6
i
I 1
hi hi
I
6I aG
I-WZL

40
 a
e
ld

erlj
0
w
z
H
I-
l C

-
F
c3
I I I I
U
4
a

oozu

41
 m
Ln
N
II
J
<
>
E
W 4-
I- I
0
Z 4
w
II
W v)
E PI
W u
I- a

W
0
Z
<
I-
cn
W
0

d
I
W
\D
6
l aJ
Ll
.-I
II
cn
a
1
s
d d 4

oozo

42
 P
I
L

W
U
Z
II

!=-I

c,
cd
v)

<
I-
m
W
a

W
h)
*
84

43
 m
m
N
I1
-I
<
>
E
w
t-
Z
H

w
t
H
I-

(D

I
W
6l
d
II
cr)
a
LI

44
 e
I
0
d

d II
m
1 P.4
W L u
td
Ip
4
m
I - - c B
d
II
v)
LL
t c,

*
m

a,
U
td
rcI
k
a,
c,
c
H

m
9)
k
1
M
d
!A

v
f 6
f b u

45
 U
I
0

I
w
hl
8
\
E (I)
4 P
11
m
a

8
4

OOZO

46
 w al
E
.rl
H &

t- [II
3
al
k
7
&
(d
k

7
al
a
E
al
f3

Il-ll-
t
dI hi dI
I
+wxa

47
rs' ai dI tsi
I 1 I

48
 ul
I
0

II

t
[II
PI
4J
(d

al
E

4J
G
H

49
 t
F
f
L

I-
I- v)
P

I I I I I I I I
t9 1
1
ni 4

uozu

50
Ln
I
W
Ea
8
4

al
$4
3
u
2
r .

al
a
E
al
E-l

51
 61
- 8
.--I

u
lu

Ln
I -d w
I \ Ll
W a,
u
c
61 .d
8
4 CJY
I! 3
0 a,
LL

Ei Ei d dI
I I I
I-WZEQ

52
 rn
PI
U
W (d
z
H
I-

53
 [I]
# PI
u
1 rd
W a
E
.rl
6
l c,
m
4
[I]
P
e
0

uozo

54
 U
(d

m' m'
I i

55
(D

I
W
m
1
4
I!
v)
a I

1111
I c
m' m'
I
m'
I
dI

56
c e r t a i n v a l u e and t h e n b e g i n s t o i n c r e a s e . Allowances f o r t h i s have

n o t been made and s i n c e t h e c o n c e n t r a t i o n of t h e l i q u i d w i l l be ever

i n c r e a s i n g , t h e m e l t i n g temperature w i l l c o n t i n u a l l y d e c r e a s e , which

can b e seen i n F i g u r e s 10, 12, 1 4 , 16, 18 and 20. It i s a l s o

i n t e r e s t i n g t o n o t e t h a t t h e v e l o c i t y of t h e i n t e r f a c e i s p r o p o r t i o n a l

t o t h e L e w i s number. The lower L e w i s numbers cause slower s o l i d i f i c a -

tion. Due t o t h e slower d i f f u s i o n r a t e s , t h e i n t e r f a c e c o n c e n t r a t i o n

i s h i g h e r and t h e r e f o r e t h e m e l t i n g temperature i s lower, r e q u i r i n g

more t i m e f o r t h e a l l o y t o r e a c h t h e m e l t i n g temperature s i n c e t h e

h e a t r a t e s remain t h e same.

F i g u r e s 21-23 i l l u s t r a t e t h e e f f e c t s of d i f f e r e n t d e n s i t i e s

i n t h e s o l i d and l i q u i d phases. The i n t e r f a c e moves a t a slower rate

when t h e d e n s i t y of t h e s o l i d i s g r e a t e r t h a n t h a t of t h e l i q u i d .

The temperature and c o n c e n t r a t i o n w i t h r e s p e c t t o t i m e does n o t change

significantly.

F i g u r e s 24-29 show t h e o v e r a l l temperature and c o n c e n t r a t i o n

distributions f o r t h e d i f f e r e n t L e w i s numbers. The c o n c e n t r a t i o n

d i s t r i b u t i o n s are r e p r e s e n t e d by t h e s o l i d l i n e s and t h e temperature

d i s t r i b u t i o n s by t h e dashed l i n e s .

I n F i g u r e 30 w e see t h e comparison of r e s u l t s o b t a i n e d from

-6
t h e a n a l y t i c a l method a t Ps = 10 ( F i g u r e 30-A) w i t h i n t u i t i v e

r e s u l t s [13] deduced from t h e phase diagram t o t h e l i m i t of z e r o

d i f f u s i v i t y i n t h e s o l i d ( F i g u r e 30-B).

57
-T I

rn
3

I
I

58
(0

I m
W P
t9
8
4

w
0
c
0
m

ai hi ai ai ai
I I- W kI
I 1 1

59
 II
t
i
l p?;
-Q
4

0
0
- .--I

I1

Q p?;
-02 !.l
0
W

0 rn
Q 3
I -io
W
W
m E
8 H
4 t-
II
cr)
a U
W
0

Q
-N

I -Q
I f
s 8
4 4

uozu

60
 G e-
0 1
0
.li
tld
a
U II

oozu 4
Ei

I-WZII
4
\
I
I n '
b I
1 I-

W
" \
\
- m
H
I3
r
H
\
\
I- \
\
\
\ 19
v) oozu

61
 h' 9 m d

6
i
I

N' Ei

N 0
6i 6i
t-wrn I-UXQ
4
I

I
\
I ' I
\
m
N
\ Ln
b
\
.--I t- 4 I-
II v, II v,
H w
w i 0 W a
r \
z
C-l
w \
t- \ t-
\ \
i \
1 19 \ t
3
J) u o z u

62
hl

-4

I-
v)
W
a

-6J
J) oozu $r
ni 6i cu' Ei

4 N
hl
8 dI 6i Ei
Ei I I
I-urn
-4

Ln
b
I-
I! v)
W
w n
r
w
I-

-61
t uozu c
ni d ni d

63
 Ln
r-
m
d Ei
I
I-wxa
-4

hl
hl
N I-
I! #
W
W 0
z
H
I-

-Ea
uozu &
d d

m d
EiI lsi
I
t-WTI1
4 -4

\
Ln I Ln
N (
b
4 I- 4 I-
II \ v) 11 v)
( H H
W \ a W a
z
H
\ r
H
I- I I-
\
\
t Ea -Ea
0 uozu t uozu
d

64
 t-wra 1
4
1
1
EJ I
Ln I I-
II I v)
I H
W 1- n
E \
H
I- \
\
\
\ \ 6l
4 uozu 4
N' d N' d

Ln hl
6l 4 4 N
d 6i
t-WZa
4
1
I
I
1
I I-
v)
\ - H
\
\
a
H \
I- \
\
\
f

Lc) u o z u 4, a

Ei

65
 I-WELL
d
1
\
m 8

m
d

It
I
I
I
I-
m
.
H
W I n
r
H
\
\
I- I
\
I j
6l
u o z u &)

(3 d

66
 A

8.5
El

I
I
I
C w '
0
N 2w
d
C B
2
H

SOLID I LIQUID
I
DISTANCE

Figure 30. Comparison of analytical method at Ps =

with intuitive results from reference [ 1 3 1 .

67
 LIST OF REFERENCES

1. American Society for Metals, Metals Handbook, Eighth Edition,

Volume 8, Metals Park, Ohio, 1973.

2. Antar , Basil N. "Solidification of a Binary Mixture", In Research

Reports--1981 NASA/ASEE Summer Faculty Fellowship Program,
G. R. Karr, J. B. Dozier, M. I. Kent and B. F. Bartfield,
Editors, NASA CR-161855, 1982.

3. Boley , Bruno A. "Time Dependent Solidification of Binary Mix-

tures". Int. J. Heat and Mass Transfer. Vol 21, No. 6,
p. 821-824, 1978.

4. Burden, R. L., Faires, J. D. and Reynolds, A. C. Numerical

Analysis. Boston: Prindle, Weber and Schmidt, 1981.

5. Carnahan, B. Applied Numerical Methods. New York: John Wiley

and Sons, 1969.

6. Carslaw, H. S. and Jaeger, J. C. Conduction of Heat in Solids.

London: Oxford University Press, 1959.

7. le Claire, A . D. "Diffusion of Metals in Metals", Progress in

Metal Physics, Chalmers, B., Editor. London: Pergamon
Press, 1961.

8. Coriell, S. R. and Cordes, M. R. and Boettinger, W. J. "Convec-

tive and Interfacial Instabilities During Unidirectional
Solidification of a Binary Alloy". Journal of Crystal
Growth. 49 (198) 1979.
9. Ince, E. L. Ordinary Differential Equations. New York: Dover,
1956.

10. Fehlberg, E. "Klassische Runge-Kutta-Formeln Funfter und

Siebenter Ordnung Mit Schrittweiten-Kontrolle." Comput-
ing. 4, p. 93-106, 1969.

11. Meyer, G. H. "On a General Theory of Characteristics and the

Method of Invariant Imbedding". SIAM J. Appl. Math. Vol.
16, No. 3, p. 488-509, May 1959.

12. . "A Numerical Method for the Two-Phase Stefan Problem".

SIAM J. Numer. Anal., Vol. 8, No. 3, p. 555-568, Sept.
1971.

13. Naumann, R. J. and Herring, H. W. Materials Processing in

Space: Early Experiments. NASA SP-443, 1980.
14. Ockendon, J. R. and Hodgkins, W. R. Moving Boundary Problems in
Heat Flow and Diffusion. London: Oxford University Press,
1975.

15. Ozisik, M. N. Heat Conduction. New York: John Wiley and Sons,
1980.

16. Reid, W. T. Riccatti Differential Equations. New York: Academic

Press, 1972.

17. Rubenstein, L. I. The Stefan Problem. Rhode Island: American

Mathematical Society, 1971.

18. Tao, L. N. "On Solidification of a Binary Alloy". Q. J. Mech.

Appl. Math. Vol. 33, pt. 2, p . 211-225, 1980.
 1. REPORT NO. 12. G O V E R N M N T ACCESSION NO. 13. R E C I P I E N T ' S CATALOG NO.
NASA CR-3962 I I
4. TITLE AND S U B T I T L E 15. REPORT D A T E
An Exact Solution for the Solidification of a Liquid March 1986
'
Slab of Binary Mixture 6 . PERFORMING ORGANIZATION CODE

7. AUTHOR(S) . 8 . PERFORMING ORGAN1 ZATION REPOR r #

B. N. Antar, F. G. Collins, and A. E. Aumalis
9. PERFORMING ORGANIZATION N A M E AND ADDRESS 10. WORK U N I T NO.
The University of Tennessee Space Institute M-513
Tullahoma, Tennessee 37388 1I . CONTRACT OR GRANT NO.
NAS8-34268
13. T Y P E OF R E P O R I & PERIOD COVEREC

I
National Aeronautics and Space Administration
Washington, DC 20546

I
15. S U P P L E M E N T A R Y N O T E S

Technical Monitor: C. F. Schafer

Prepared for Systems Dynamics Laboratory, Geo. C. Marshall Space Flight Center
16, ABSTRACT

The time dependent temperature and concentration profiles of a one dimensional

finite slab of a binary liquid alloy is investigated during solidification. The
governing equations are reduced to a set of coupled, nonlinear initial value
problems using the method outlined by Meyer in ref [11,12]. Two methods will be
used to solve these equations. The first method uses a Runge-Kutta-Fehlberg
integrator to solve the equations numerically. The second method comprises of
finding closed form solutions of the equations.

17. KE? WORDS 18. D l S T R I B U T I ON S T A T E M E N T

Correction Subject Category: 44
Binary Mixture
Solidification Unclassified-Unlimited
Fluid Mechanics
Stefan Problem

19. S E C U R I T Y C L A S S I F . (ofthln ropart) 20. S E C U R l T Y C L A S S I F . (of thln page) 21. NO. O F PAGES 22. PRICE

Unclassified Unclassified 76 A05

NASA-Langley, 19&