...

mbp01.pdf

by user

on
Category: Documents
2

views

Report

Comments

Description

Transcript

mbp01.pdf
Moving Boundary Problems
Ernesto Gutierrez-Miravete
Rensselaer at Hartford
April 28, 2006
1 Introduction
A substantial proportion of problems involving heat, mass or momentum transfer commonly
encountered in engineering practice require nding the solution inside a domain whose boundary is unknown and may even be a function of time. Free boundary problems are obtained
under steady state conditions whereas transient situations give rise to moving boundary
problems. This presentation reviews the general formulation of moving boundary problems
in conduction heat transfer with change of phase and in solid state diusion with phase transformations. The few existing analytical solutions are rst introduced and examined using the
symbolic manipulation code Maple. Successful and robust numerical solution methodologies
are then presented and demonstrated using simple Fortran programs. Finally, a sample of
the industrial application of the above ideas, drawn from the author's experience, will be
presented and discussed.
2 Moving Boundary Problems in Heat Conduction
Changes in atomic arrangements in condensed phases involve the release or absorption of
thermal energy at the phase interface. Common examples are the release of latent heat of
freezing during the solidication of a metal, the austenitization of steel during heat treatment
and crack sealing during wide gap braze repair of airfoils. The mathematical formulation of
these problems require the statement of heat equations for the phases involved, specication
of necessary initial and boundary conditions and a description of conditions at the interface
between the phases. Since the position of the interface is not known in advance and must
be determined as part of the solution, the problem is nonlinear.
2.1 Mathematical Formulation of Freezing and Melting Problems
Consider the problem of solidication in a semiinnite region 0 < x < 1 initially liquid
and at uniform temperature T (x; 0) = Ti > Tm where Tm is the melting/freezing point of
1
the substance in question and subjected to a xed temperature T (0; t) = T0 < Tm at its
boundary x = 0. After some time t, a solid skin has formed and the interface separating
solid and liquid phases (t), moves along the positive x direction.
The formulation of the problem requires nding functions Ts(x; t), Tl (x; t) and (t) such
that
@ 2 Ts 1 @Ts
=
in 0 < x < (t),
in (t) < x < 1, subject to
@x2
s @t
@ 2 Tl
@x2
l
= 1 @T
@t
l
Ts (0; t) = T0
at x = 0,
Tl (x; t) ! Ti
as x ! 1, and
Ts (; t) = Tl (; t) = Tm
and
@T
d
@Ts
kl l = L = Lv
@x
@x
dt
x = (t). The rst of the last
ks
at the interface boundary
two conditions species perfect
contact at the solid-liquid interface while the second is a statement of the thermal energy
balance there. Here L is the latent heat of freezing (J/kg) and v = d=dt is the velocity of
advancement of the solid-liquid interface.
If the densities of solid and liquid phases are dierent, (say s > l) the s must be
substituted for in the RHS of the interface heat balance equation. Also, if convection in
the liquid is the dominant mode of heat transfer near the interface, the heat ux in the liquid
must be given instead by h(T1 Tm).
In 3D systems, the corresponding form of the heat balance at the interface is
ks
@Ts
@n
kl
@Tl
@n
= Lvn
where vn is the interface velocity in the normal direction.
2
2.2 Exact Solutions
Consider the problem of melting of a semi-innite (0 < x < 1) solid, which is initially at a
uniform temperature Ti Tm . Suddenly, the temperature at x = 0 is xed at T0 > Tm and
maintained there. The solid begins to melt and the liquid-solid interface advances along the
positive x direction. since Ts(x; t) = Tm one must nd Tl (x; t) in 0 < x < (t) such that
@ 2Tl 1 @Tl
=
@x2
subject to
l @t
Tl (0; t) = T0
Ts (x; 0) = Ti
for x > 0 and t = 0. The conditions at the interface x = (t) are
Tl (; t) = Tm
and
d
@T
kl l = L
@x
dt
A solution of the form
x
Tl (x; t) = T0 + B erf [ p ]
2 t
l
where the constant B is still to be determined satises the heat equation and the boundary
condition at x = 0. Satisfaction of the perfect thermal contact condition at the liquid-solid
interface requires that
B=
Tm T0
p
erf [=2 l t]
m T0
= Terf
()
where = =2pl t must be determined so as to satisfy the remaining heat balance condition
at the interface. The required solution thus becomes
p
Tl (x; t) T0 erf [x=2 l t]
= erf ()
Tm T0
Finally, the required value of is the root of the trascendental equation
C (T T )
2
e erf () = p 0p m
L 3
A related problem is the solidication of a large mass of supercooled liquid (i.e. Ti < Tm)
from its surface maintained at Tm. Furthermore, no heat loss is allowed through the solid so
that its temperature is everywhere equal to Tm.
By an argument similar to the one above, the temperature prole in the supercooled
liquid is given by
p
Tl (x; t) Ti erfc[x=2 l t]
= erfc()
Tm Ti
where the required value of is the root of the trascendental equation
C (T
T)
2
e erfc() = p mp i
L Consider again the initial problem of a semiinnite liquid (0 x < 1), initially at
uniform temperature T (x; 0) = Ti > Tm whose surface at x = 0 is maintained at T (0; t) =
T0 < Tm begining at t = 0. A solid shell will form and the solid-liquid interface will advance
along the positive x direction. Since here T0 < Tm and Ti > Tm temperature proles will
develop in the two phases. Solutions for Ts(x; t) and Tl (x; t) of the forms
p
Ts (x; t) = T0 + A erf [x=2 s t]
and
p
Tl (x; t) = Ti + B erfc[x=2 l t]
were A and B are constants yet to be determined, satisfy the heat equation, the boundary
conditions at x = 0 and x ! 1 as well as the initial condition.
The condition of perfect contact at the solid-liquid interface requires that
A=
and
B=
Tm T0
erf ()
Tm q Ti
erfc[ s =l ]
where
=
(t)
p
2 t
s
4
with the above, the required solutions are
p
Ts (x; t) T0 erf [x=2 s t]
= erf ()
Tm T0
and
p
Tl (x; t) Ti erfc[x=q2 l t]
=
Tm Ti
erfc[ s =l ]
Finally, a trascendental equation for the determination of is obtained by forcing the
above expressions to satisfy the interface heat balance, the result is
p
( = )
kl s 1=2 Tm Ti e 2q
L e 2
+ ( )
=
erf () ks l Tm T0 erfc[ s =l ] Cp;s (Tm T0 )
s
l
Finally consider the problem involving a semiinnite medium (0 < x < 1), initially
liquid at constant temperature Ti > Tm. The properties of the liquid are kl , l , Cp;l and
those of the associated solid ks, s, Cp;s. This medium is then brought at time t = 0 into
perfect thermal contact with another semiinnite medium ( 1 < x < 0) with properties k0,
0, Cp;0 and initially at constant temperature T0 (x; 0) = 0 < Tm . This may be a model of the
initial stages of solidication of a casting in a mold. Temperature distributions in the mold,
solidied shell and in the liquid develop. As before, the problem can be solved assuming a
specic form for the temperature distributions and then determining the constants involved
by enforcing the initial and boundary conditions. The nal result is
T0 (x; t) =
for 1 < x < 0
Ts (x; t) =
ks 01=2Tm
(1 + erf [ 2px t ])
1
=2
1
=2
0
ks 0 + k0 s erf ()
Tm
(ks01=2 + k0s1=2erf [ 2px t ])
1
=2
1
=2
s
ks 0 + k0 s erf ()
for 0 < x < (t) (the solidied shell), and
(Ti Tm) erfc[ px ])
Tl (x; t) = Ti
erfc[(s =l )1=2]
2 l t
p
where = =2 l t is now the root of
k0 s1=2e 2
kl s1=2 (Ti Tm )e 2 =
L1=2
=
ks 01=2 + k0 s1=2erf () ks l1=2 Tm erfc[(s =l )1=2] Cp;s Tm
s
5
l
2.3 Numerical Methods
Numerical solution techniques for phase change problems must account for the changing
location of the phase boundary. Two commonly used procedures are variable grid methods
and xed grid techniques. In variable grid methods one makes the (moving) phase boundary
coincide with a particular mesh point at all times, while in xed grid techniques the boundary
location is obtained by interpolation from the computed values of temperature on a xed
mesh.
Since xed grid methods are usually easier to implement we shall focus on those here. The
most commonly used xed grid method for the numerical solution of solidication problems
is based on the enthalphy-temperature formulation.
Consider a slab of liquid material, initially at T (x; 0) = Ti > Tm which is made to solidify
by xed the temperature at its boundaries T (0; t) = T (l; t) = T0 < Tm. Instead of writing
separate equations for the solid and liquid phases, the heat equation is written as follows:
@H
@t
@ @T
= @x
(k @x )
where the enthalpy H (J/kg) is a step function of temperature, i.e.
(
H = sCCp;s((TT TTm));+ TL;< TTm > T
l p;l
m
s
m
Introduce now a xed mesh in space (mesh spacing x) with N + 1 nodes located at
xi where i = 1; 2; :::; N + 1 and a mesh in time (mesh spacing t) with nodes tj where
j = 1; 2; 3; :::.
Discretization of the heat equation by means of a simple explicit nite dierence or nite
volume scheme yields
Hi;j +1
Hi;j
ki+ 21 ;j T +1x T
i
;j
i;j
ki
x
1 ;j
2
Ti;j Ti 1;j
x
=
t
where the two subscripts on the variables are needed to specify the mesh point where the
quantity is to be evaluated.
The numerical calculation is carried out from an initial thermal condition by rst computing values of enthalpy for all spatial nodes at the new time step (Hi;j+1) explicitly from
the last expression. Temperature values for all spatial nodes at the new time step are then
calculated from the given enthalpy-temperature relationship. The position of the moving
boundary is determined by locating spatial nodes around Tm and interpolation. The calculation is then repeated to advance to the next time step. Since this is an explicit scheme,
the mesh spacings must fulll the CFL condition, namely
2
t = 2x
6
where is the smaller of s and l .
Alternatively, an implicit scheme may be employed but this will require implementation
of an iterative process to advance the solution each time step.
2.4 Exercises
Exercise 1. The surface of a large solid Titanium mass initially at Ti Tm = 1933K ,
is suddenly raised to a temperature of T0 = 2000K and maintained there. Estimate the
advancement of the melting interface. Consider k = 15:5W=mK , = 4; 505kg=m3, Cp =
522J=kgK , L = 392; 648J=kg.
Exercise 2.
The surface of a large supercooled liquid Nickel mass initially at Ti = 1528K ,
is suddenly seeded with a crystal nucleus which induces solidication of the melt. Estimate
the advancement of the solidication interface. Consider k = 69:1W=mK , = 8; 847kg=m3,
Cp = 418J=kgK , L = 293; 065J=kg, Tm = 1728K .
Exercise 3.
Use the properties of the error function to verify that the given expression
for Ts in the problem of solidication of a semiinnite melt satises the heat equation as
well as the boundary condition at x = 0 while that for Tl satises the heat equation, the
boundary condition at x ! 1 as well as the initial condition.
Exercise 4.
The surface of a large molten mass of Aluminum initially at Ti = 1033K ,
is suddenly lowered to T0 = 333K and maintained there. Estimate the advancement of
the solidication interface. Consider k = 238W=mK , = 2; 700kg=m3, Cp = 945J=kgK ,
L = 389; 177J=kg, Tm = 933K .
Exercise 5.
The surface of a large molten mass of Aluminum initially at Ti = 1033K ,
is suddenly placed in intimate thermal contact with a massive steel chill originally at T0 =
333K . Estimate the advancement of the solidication interface. Consider the following
property values for steel k = 52:3W=mK , = 7; 900kg=m3, Cp = 470J=kgK .
Exercise 6.
The temperature of the surface of a large molten mass of Iron initially at Ti =
1910K , is suddenly lowered to T0 = 25K and maintained there. Estimate the advancement of
the solidication interface using an explicit nite dierence method. Consider k = 30W=mK ,
= 7; 300kg=m3, Cp = 600J=kgK , L = 277; 400J=kg, Tm = 1810K .
7
3 Moving Boundary Problems in Solid State Diusion
At equilibrium, the values of chemical potential of each substance in a multicomponent
system are the same in all constituent phases of the system. For instance, in a binary system
consisting of two substances A and B and existing in two phases and , the equilibrium
conditions are
A = A
and
B = B
Sometimes it is more convenient to state the equilibrium conditions in terms of the
activities of the various substances, i.e.
aA = aA = aA
and
aB = aB = aB
If for some reason, gradients in the chemical potentials of constituent substances develop,
the system naturally changes in such a way that the gradients disappear. The mechanism by
which chemical potential gradients are reduced is atomic transport of substance from regions
of high potential to regions of low potential. This is entirely analogous to the situation
encountered in heat conduction where thermal energy is transported from regions at high
temperature toward regions at low temperature in a natural attempt by the system to reduce
temperature gradients. .
An interesting and commonly encountred situation in multicomponent multiphase systems in equilibrium is the coexistence of dierent phases with dierent concentrations of
substance (at the same chemical potential). For instance, for the same binary system above
the equilibrium concentrations of the two components may be given by
cA = cA 6= cA = cB
and
cB = cB 6= cB = cB
Therefore, if such as system is taken out of equilibrium thus introducing gradients in
the chemical potentials of the constituent substances, the natural response of the system is
the atomic displacement (diusion) of the constituent substances in an attempt to bring the
system back into equilibrium conditions. In this process the relative amounts of phases and may change with time. This is a phase transformation process and a typical moving
boundary problem, the boundary being the interface separating the two coexisting phases.
8
3.1 Mathematical Formulation of Moving Boundary Problems in
Solid State Diusion
The mathematical formulation of these problems require the statement of diusion equations
for the phases involved, specication of necessary initial and boundary conditions and a
description of conditions at the interface between the phases. Since the position of the
interface is not known in advance and must be determined as part of the solution, the
problem is nonlinear.
Consider the problem of phase transformation in a semiinnite region 0 < x < 1 of
a binary system initially all and at uniform potential of the constituent susbtances. We
will consider only the case where the atomic mobilities of the component substances are
signicantly dierent so that only the diusion of the fastest moving species needs to be
examined and the subscripts A and B can be omitted. Therefore, let the activity and
concentration of this substance be a and c, respectively. Let the initial condition of the
system be given by a (x; 0) = ai > a where a is the value of the activity of the substance
in question at equilibrium of the phases and (in terms of concentration c (x; 0) = ci).
We shall also assume that although a = a = a, the concentrations c = c < c = c.
Now let the activity at its boundary x = 0 be xed to a value a0(0; t) < a (concentration
c(0; t) = c0 < c ). As a result, a layer of phase forms at the surface and the interface
separating the and phases (t), moves along the positive x direction. It is commonly
assumed that equilibrium conditions are quickly reached at the interface so that there (and
only there) c = c and c+ = c. This is the condition of local equilibrium at the interface.
The formulation of the problem requires nding functions c(x; t), c (x; t) and (t) such
that
@ 2 c
= 1 @c
in 0 < x < (t),
in (t) < x < 1, subject to
at x = 0,
@x2
D @t
@ 2 c
@x2
= D1 @[email protected]
c (0; t) = c0
c ! ci
as x ! 1.
At the interface x = (t) one has the local equilibrium condition
a (; t) = a (; t) = a
9
(or, equivalently
and the interfacial mass balance
D
@c
@x
D
c = c; c = c
@c
@x
= (c
c)
d
dt
= (c
c)v
Here (c c) is the jump in concentration at the interface and v = d=dt is the velocity
of advancement of the interface.
In 3D systems, the corresponding form of the mass balance at the interface is
@c
@n
@c
@n
= (c c)vn
where vn is the interface velocity in the normal direction.
D
D
3.2 Exact Solutions
Consider the situation involving a semi-innite sample of a binary system in single phase ()
state with a uniform initial concentration of diusing substance c(x; 0) = ci. At t = 0 the
surface concentration is increased and therein maintained at a xed value c(0; t) = c0 > c
(surface enrichment). As a result, a layer of phase forms at the surface and its thickness
increases with time with the interface moving along the positive x direction.
The mathematical formulation of the problem requires nding functions c(x; t), c (x; t)
and (t) such that
@ 2 c
= 1 @c
in 0 < x < (t),
in (t) < x < 1, subject to
at x = 0,
as x ! 1, and
@x2
D @t
@ 2 c
@x2
= D1 @[email protected]
c (0; t) = c0
c (x; t) ! ci
a (; t) = a (; t) = a
10
(or, equivalently
c = c; c = c
and
@c
@c
D
@x
@x
at the interface boundary x = (t).
D
= (c
c )
d
dt
This problem is readily solved by assuming the solution can be expressed in terms of
error functions. A solution of the form
x
c (x; t) = c0 + B erf [ p ]
2 Dt
where the constant B is to be determined, satises the diusion equation inside the layer
as well as the boundary condition at x = 0.
Satisfaction of the equilibrium condition at the interface requires that
B=
c0 c
erf ()
p
where = =2 D t must be determined so as to satisfy the remaining condition at the
interface.
Inside the phase the appropriate form of the solution is
x
c (x; t) = ci + C erfc[ p ]
2 Dt
where the constant C is determined by invoking the equilibrium condition at the interface.
Finally the required solution for the entire system thus becomes
p
c (x; t) c0 erf [x=2 D t]
= erf ()
c c0
in the layer and
p
c (x; t) ci
erfc[x=2 D t]
= erfc((D =D)1=2
c ci
in the interior.
The interface position is given by
p
(t) = 2 D t
11
Finally, the required value of is the root of the trascendental equation
c0 c
(
c ci )e 2
2
p
p 1=2erfc( 1=2) = c c
e
erf ()
where
= DD
The above solution can be applied to estimate the progress of oxidation of a pure metal
(ci = 0) where is the oxide layer and oxygen is the diusing substance.
Finally consider the problem involving a semiinnite medium (0 < x < 1), initially all
consisting of a pure substance B (i.e. cB = c = 1 in atomic fraction units). This medium
is then joined at time t = 0 with another semiinnite medium ( 1 < x < 0) consisting of
pure A (i.e. cA = 1; cB = c = 0) in phase . The two phases coexist in equilibrium over
the concentration range from c to c. Finally, assume that the diusivity
of substance B
is much larger than that of A in both phases so that DB = D and DB = D .
This is a model of a diusion couple, an experimental system commonly used to investigate diusional processes in solids. The problem can be solved as before by assuming a
specic (error function) form for the concentration distributions and then determining the
constants involved by enforcing the initial and boundary conditions. The nal result for
1 < x < 0 is
x
c
[1
+
erf ( p )]
c (x; t) =
1 + erf ()
2 Dt
and for 0 < x < 1,
1
erf (= 1=2 ) + (1 c )erf ( px )]
c (x; t) =
[
c
erfc(= 1=2 )
2 D t
Further, the interface position is given by
p
(t) = 2 D t
where is the root of
ce 2
(1 cp) 1=2e 2=
p
(c c) (1 + erf ()) (c c) erfc(= 1=2) = 1
where
= DD
12
3.3 Numerical Methods
Numerical solution techniques for diusion problems with phase change must account for
the changing location of the phase boundary. As in the case of heat transfer, variable and
xed grid methods can be used. Since xed grid methods are usually easier to implement
we shall again focus on those here. The xed grid method used for the numerical solution
of moving boundary problems is based on the activity formulation.
Consider a slab of binary material in phase with initial concentration of diusing
substance c(x; 0) = ci > c. At time t = 0, the concentration of diusing substance at the
boundaries of the slab is lowered and maintained at c(0; t) = cS < c. As a result, layers of
form at the surface and interfaces move towards the center of the slab. Instead of
writing separate equations for the solid and liquid phases, the diusion equation is written
as follows:
@c
@t
@ @a
= @x
(D @x )
where the concentration c is a step function of activity, i.e.
( + c ; a < a
c = ((aa aa )) +
c; a > a
where c = c c.
Introduce now a xed mesh in space (mesh spacing x) with N + 1 nodes located at
xi where i = 1; 2; :::; N + 1 and a mesh in time (mesh spacing t) with nodes tj where
j = 1; 2; 3; :::.
Discretization of the heat equation by means of a simple explicit nite dierence or nite
volume schemes yields
ci;j +1
ci;j
Di+ 21 ;j a +1x a
i
;j
i;j
Di
x
1 ;j
2
ai;j ai 1;j
x
t =
where the two subscripts on the variables are needed to specify the mesh point where the
quantity is to be evaluated.
The numerical calculation is carried out from an initial concentration condition by rst
computing values of concentration for all spatial nodes at the new time step (ci;j+1) explicitly
from the above expression. Activity values for all spatial nodes at the new time step are then
calculated from the given concentration-activity relationship. The position of the moving
boundary is determined by locating spatial nodes around a and interpolation. The calculation is then repeated to advance to the next time step. Since this is an explicit scheme, the
mesh spacings must fulll the CFL condition, namely
2
t = 2Dx
13
where D is the smaller of D and D .
Alternatively, an implicit scheme may be employed but this will require implementation
of an iterative process to advance the solution each time step.
3.4 Exercises
Exercise 1. Use the properties of the error function to verify that for the two phase
binary system , the given expression for c satises the diusion equation as well as the
boundary condition at x = 0 while that for c satises the diusion equation, the boundary
condition at x ! 1 as well as the initial condition.
Exercise 2. The surface of a large piece of pure iron is exposed to a carburizing atmosphere which xes the surface concentration of carbon to c0 = 1 percent. At the treatment
temperature of 900 degrees Celsius a carbon rich austenite ( ) phase forms on the surface of
the carbon lean ferrite () core. Estimate the advancement of the interface. Consider
D = 5 10 12, D = 5 10 10 (both in m2 =s), c = 0:001, c = 0:01 (both in weight
percent).
Exercise 3. A large sample of iron-carbon alloy contains 1 percent carbon. The sample is
heated to 900 degrees Celsius (whence it becomes single phase - austenite or ) and exposed
to an environment that lowers the surface concentration of carbon to 0. As a result, a layer
of phase forms on the surface and the interface moves towards the interior of the
sample. Estimate the advancement of the interface using an explicit nite dierence method.
4 References
1.- L.I. Rubinstein, The Stefan Problem, American Mathematical Society, Providence,
1971.
2.- J.R. Ockendon and W.R. Hodgkins, Moving Boundary Problems in Heat Flow and
Diusion, Clarendon Press, Oxford, 1975.
3.- J. Crank, Free and Moving Boundary Problems, Clarendon Press, Oxford, 1984.
4.- V. Alexiades and A.D. Solomon, Mathematical Modeling of Melting and Freezing
Processes, Hemisphere, Washington, 1993.
5.- J. Philibert, Atom Movements: Diusion and Mass Transport in Solids, Les Editions
de Physique, France, 1991.
6.- R. Sekerka, C.L. Jeanls and R.W. Heckel, The Moving Boundary Problem, Chapter
IV in H. I. Aaronson (ed), Lectures on the Theory of Phase Transformations, The Metallurgical Society, New York, 1975.
14
Fly UP