Crank J. Free and Moving Boundary Problems
Подождите немного. Документ загружается.
112
Analytical solutions
(1929) solved
the
problem
of
solidification in x
>0
with initial tempera-
ture
U above
the
melting temperature
Urn
and
the
surface x = 0 main-
tained
at
zero for t >
O.
Carslaw
and
Jaeger (1959, p. 293) give
the
following alternative presentation of Lightfoot's method.
The
tempera-
ture
due
to
the
plane moving
heat
source (3.24) may
be
written
u(x
t) =
L,Lt
S(7')
d-r
[e-{x-s(T)}2/4k(t-T)
-
e-{x+s(T)}2/4k(t-T)J,
(3.25)
,
2c('1Tk)'i
(t-T)~
which
is
obtained by summing contributions
to
the
temperature
at
(x, t)
from sources of strength LpS(T) dT instantaneously
at
x =
S(T)
for
the
short time interval T to T +
dT,
by integrating over
the
whole period
of
time T=O
to
t.
In
addition,
the
temperature
at
(x, t) is contributed
to
by
the
heat
originally in
the
region x > 0 associated with
the
initial uniform tempera-
ture
U.
This contribution, having in mind
the
surface condition of zero
temperature
at
x = 0, is
x
w(x, t) = U
erf
2(kt)!
'
(3.26)
which also satisfies
the
initial condition w(x, 0) =
U,
x>O. Thus,
the
condition
that
the
temperature
at
x = s(t) is
the
melting temperature
Urn
is,
on
combining (3.25)
and
(3.26),
u{s(t),
t}+
w(s(t), t) =
Urn,
(3.27)
which is
an
integral equation for s(t).
Lightfoot solved
the
problem by adopting
the
similarity transformation
s(t) = 2A(kt)f (3.28)
and
expressing
the
integral (3.25) in terms
of
error
functions.
In
terms
of
new variables
(3.25) becomes
y = x/s(t),
=
(1_Z2)2
T 1
2'
+Z
2ll
u(x, t)
=-1
(1
1
-
Iz),
C'1T
2
where
the
integrals /1
and
12
reduce
to
/1
=l'1Te
A
'
erfc
A(1
+erf
AY),
=l'1Te
A
'
erfc Ay(1
+erf
A),
12
=
l'1Te
A
'
erfc
AY
erfc
A,
y<1,
y>1,
(3.29)
(3.30)
(3.31)
(3.32)
(3.33)
Similarity solutions and moving heat sources
113
where clearly from
the
first of (3.29), y
~
1 correspond
to
the
liquid and
solid regions respectively,
and
y = 1 is
the
solidification front where (3.31)
and
(3.32)
are
identical. Substituting (3.31)
and
(3.33) with y = 1, and
(3.26),
the
integral equation (3.27) becomes
Lhr!
--
e
A2
erfc A erf A + U erf A =
Urn,
c
(3.34)
which is
the
Neumann expression for A in this problem when
Kl
= K
z
and
kl = k
z
.
Lightfoot (1929) extends (3.25) so as
to
include an infinite
number of image terms instead of just the two given and obtains an
approximate solution for
the
solidification of a finite plane sheet.
Rathjen and Jiji (1971) extend Lightfoot's idea of treating
the
solidifi-
cation front as a moving heat source
to
obtain an analytical solution of
the
two-dimensional problem of
the
solidification of a liquid, initially
at
a
uniform temperature
and
occupying
the
quarter-space
x,
y
>0,
subject to
a constant temperature on
the
boundaries x = 0, y =
O.
They express
the
temperature as
the
sum of two parts.
The
first is a standard solution of
the
heat
flow equation in a rectangular comer, without phase changes,
satisfying
the
initial
and
boundary conditions of
the
present problem;
the
second
part
represents
the
contribution from
the
moving heat source by
using Green's functions for
the
rectangular
comer
with zero initial and
surface temperatures.
After
some manipulation and use of similarity
variables
xlt!,
ylt~,
an
integro-differential equation for
the
movement of
the
solidification boundary is obtained. Numerical solutions are computed
by using superhyperbolae
to
approximate
the
position of
the
moving
boundary. Values of
the
temperature field follow for all times. Compari-
sons
are
made with finite-difference calculations and with experimental
results.
By a convenient choice of temperature variables, Lightfoot's restriction
of equal conductivities and specific heats in
the
two phases
is
replaced by
the
requirement
that
only
the
ratios of conductivity to specific heat need
be
the
same. Even this restriction
is
unnecessary if
the
liquid is initially
at
the
solidification temperature.
Budhia
and
Kreith (1973) proceed very much in
the
same way as
Rathjen
and
Jiji (1971)
to
obtain analytical solutions for melting
or
freezing of materials in a wedge-shaped space with wedge angle between
o and 360 degrees.
The
material is initially
at
a uniform temperature
within
the
wedge
and
the
surface temperatures are kept constant
but
not
necessarily equal.
The
solutions are expressed in radial coordinates
(r,
«»
and time
t,
and
a similarity transformation
rlt!
is
justified.
To
solve their
integro-differential equation Budhia
and
Kreith represent
the
moving
boundary by a one-parameter hyperbola.
Computed results
are
given for three examples:
(i)
liquid
at
the
fusion
114
Analytical solutions
temperature inside a 90° wedge with
the
surface temperatures equal
and
below
the
fusion temperature;
(ii)
the
same
problem
except
that
the
liquid
is initially
at
a temperature higher than
the
fusion temperature
and
wedge
angles
of
60°
and
270° are considered as well as 90° (a way
of
improving a
first analytic approximation
is
used when
the
diffusivities in
the
two
phases differ);
and
(iii)
the
same problem
but
the
surfaces
of
the
90°
wedge
are
maintained
at
unequal temperatures lower
than
the
fusion
temperature.
3.5. Integral-equation
formulations
The
use
of
moving
heat
sources in §3.4 is a special way of deriving
an
integral equation
and
likewise
the
similarity solution is a particular form.
Several
other
more
general ways
of
formulating moving-boundary
prob-
lems in terms
of
integral equations
and
their solutions have
been
found
useful.
One
advantage
of
introducing integral equations in linear heat-
flow problems generally
is
that
only
the
values
of
the
unknowns
on
the
boundaries
of
the
domain
enter
into
the
formulation. Although Stefan-
type problems
are
essentially non-linear because
of
the
condition
on
the
moving boundary, nevertheless linearity frequently exists in each
of
the
domains
on
either
side
of
the
boundary
ano
the
integral equation method
is still useful.
3.5.1. Integral transforms
Evans
et
al. (1950) used Laplace transforms
to
express a one-phase
moving boundary
problem
of
the
kind described in §1.2.1
but
with a
constant
heat
source applied
to
the
outer
surface
of
the
domain. In
the
nomenclature
of
Evans
et
al. (1950)
the
problem is
to
find x = x(t) where
a
2
=K/cp,
u
=0,
x
~x(t),
A
d~~t)
=
::
(x(t), t), A =
Lp/
K,
x(O)=O,
au/ax
= g(t),
x=O,
t~O,
(3.35)
(3.36)
(3.37)
(3.38)
(3.39)
where g(t)
is
a time-dependent
heat
input
and
L
here
is
the
latent
heat
of
recrystallization from a to
(3
crystals when
iron
is
heated
through this
phase-change temperature.
In
order
to obtain
the
integral equation satisfied by
the
moving bound-
ary, x = x(t), we rewrite this relationship as t =
f(x)
and
reformulate
the
Integral-equation formulations
115
problem as
au 2 a
2
u
at
=a
ax2'
t>f(x)
u(x,t)=o,
t~f(x),
au/ax = g(t),
x=O,
t;;;:oO,
f(O)
=0,
au A
ax (x, f(x)) =
f(x)
,
f=df/dx.
Define
the
Laplace transform
of
u(x, t) as
(3.40)
(3.41)
(3.42)
(3.43)
(3.44)
u(x,p)=
foou(x,t)e-ptdt=j"" ue-Ptdt, (3.45)
t
f(x)
and
hence from (3.40) obtain
d
2
u _ pu =
-Ae-pf(x)
dx
2
a
2
'
with
the
boundary conditions
au/ax=g,
x=O;
l~
u(x,
p)
=
O.
The
solution
-of
(3.46) is
(3.46)
(3.47)
u(x, p) =
~~
LOOex
p
{
-Pf(~)}[exp{
-~~
Ix-~I}+exp{
-~!
(x+~)
}]d~
ag
{_pi}
--exp
-x,
p!
a
(3.48)
and
the
inverse transform is
(3.49)
116
Analytical solutions
Use
of
the
condition
u(x,
f(x))
= 0 from (3.41)
in
3.49) gives
i
f(X)
g(t-T)
(_X2)
2
--exp
--
dT
Tf
4a
2
T
,
i""
1
[{
_(X_~)2
}
=A
{f(x)-f(~)}f
exp
4a2[f(x)-f(~)]
+
{
-(x+~?
}]
+exp
4a
2
[f(x)-fW]
d~,
f(x)
>
f(~),
(3.50)
which is
an
integral equation for
the
determination
of
t =
f(x).
In
order
to
solve equation (3.50) we
return
to
the
specification of
the
moving
boundary as x =
x(t)
and
put
t =
f(x),
T =
f(~),
x =
x(t),
~
=
~(T)
which is
permissible provided t =
f(x)
, is monotonically increasing.
After
these
changes
on
the
right-hand side
of
(3.50) we replace
the
new
t-T
by
T'
and
finally drop
the
superscript
and
write
T'
= T so
that
(3.50) becomes
2
it
g(t-:-T)
exp{-X
2
(t)} dT
T2
4a
2
T
=A
it
X(t:T)
[exp{-[X(t)-x(t-T)f}
+
T~
4a
2
T
{
-[X(t)
+
X(t-T)]2]
+exp
4 2
dT.
aT
T>O. (3.51)
For
the
case
of
g(t) = g = constant, Evans et
al.
(1950) obtain
the
first two
terms in a power series solution for x =
x(t)
i.e.
(3.52)
where
(3.53)
which agrees with
the
first two terms of a Taylor series obtained by
the
same authors for
x(t)
about t =
O.
When
g(t)
is
an
analytic function of
t,
the
first two coefficients
Cl>
C
2
in
the
solution
of
(3.51)
are
also shown
to
agree with
the
corresponding terms of
the
appropriate Taylor series. They
are
C
1
= g(O)/A, (3.54)
Ockendon (1975) shows how
an
integral formulation of
the
oxygen
diffusion problem described by equations (1.63) in §1.3.10 follows by a
similar use
of
Laplace transforms.
The
integral equations
are
found
to
be
those derived
by
Hansen
and
Hougaard (1974)
and
they
are
obtained in
Integral-equation formulations
117
§3.5.2 by using Green's functions. Again it
is
necessary
that
x(t)
should
be
monotonic.
Ockendon (1975) also draws attention to the usefulness of Fourier
transforms in infinite
or
semi-infinite domains and cites
the
Fourier
transform derivation of Lightfoot's (1929) integral equation based
on
a
moving
heat
source as in equations (3.25), (3.26), and (3.27).
3.5.2. Green's /unctions
The use
of
Green's functions for the solution of heat flow problems
subject to conditions prescribed on fixed boundaries
is
well known.
Carslaw and Jaeger (1959, Chapter 14) give the basic theory and solu-
tions
to
a selection of standard problems in terms of appropriate Green's
functions. Kolodner (1956, 1957) formulated Stefan problems in terms of
integral equations, including the problem of the freezing of a lake of finite
depth, by using simple Green's functions, though
he
does not expljcitly
say so.
It
is
interesting
that
in such an early paper Kolodner remarked
that
the
equations in cases of physical interest are all of Volterra
type
of
the
second kind and hence are amenable to numerical treatment. Rubin-
stein (1971) based an analysis of the existence, uniqueness, and stability of
solutions of Stefan problems on integral equation formulations incor-
porating Green's functions.
He
also briefly examined the numerical
evaluation of integral solutions coupled with estimates of accuracy. More
recently Rubinstein
(1980a) discussed the application of integral equation
techniques to several Stafan problems including a one-din'lensional prob-
lem with analytical input data,
the
solidification of a binary alloy, an
axially symmetric problem with concentrated thermal capacity (see
§1.3.11), and a problem with a hyperbolic instead of the more usual
parabolic
heat
conduction equation. Collatz (1978) too used integral
equations
tc?
find practical error bounds. Chuang and Szekely (1971)
employed Green's functions to solve, in integral equation form, the
problem of a solid slab, symmetrically placed in its own melt. There
is
a
prescribed initial temperature, not necessarily uniform, and
on
the mov-
ing,
outer
surfaces x =
±X(t)
the temperatures,
U,
are the melting temper-
ature,
UM,
and there
is
a second condition of the form
au
a-+
(3u
=
t(t),
ax
x =
±X(t)
,
with
X(t)
= e at t =
O.
The
numerical evaluation of the integral solution
proceeds by a process of successive approximation in which the locus of
the melting boundary
is
approximated by linear segments in the (x, t)
plane.
Later
papers by Chuang and Szekely (1972) and Chuang and
Ehrich (1974) deal with corresponding cylindrical and spherical problems.
118
Analytical solutions
An
extension to the solidification of iron-carbon alloys
is
outlined by
Chuang et
al.
(1975), who give some numerical results.
Hansen and Hougaard (1974) applied Green's functions to obtain
solutions of the oxygen diffusion problem prescribed by equations (1.57)
to (1.60) in §1.3.10. This paper contains useful details of the analysis and
also numerical comparisons with those obtained
by
Crank and
Gupta
(1972a,b).
It
has
the
added interest that Hansen and Hougaard treat
the
problem in its implicit form (Sackett 1971) in which
the
moving boundary
condition does not contain the velocity of the boundary explicitly. Their
solution proceeds along the following lines:
A Green's function G defined by
1
{(-(x
-
X')2)
(-(X
+
X')2)}
G(x, x', t' -
t)
exp + xp
--'-,..----c"-
27T.J
7T(t'
-
t)
4(t' -
t)
e 4(t' -
t)
(3.55)
for
t'>t
and
G(x,x',
t'-t)=O
(3.56)
for
t > t', satisfies
;PG
aG
( ')
('
)
-+-=-8
x-x
8
t-t
ax
2
at
'
(3.57)
where 8
is
the
Dirac delta function, for x and x'
~O
and for all
t.
It
also
satisfies the boundary condition
G,,(O,
x',
t'-
t)
= 0 (3.58)
for all
t and x'>O.
Following the standard treatment of Green's functions the integral
I
t'
l"o(t)
[
~2C
ac)
(a
2
G
a~]
G
---
-c
-+-
dxdt
o x
2
at
ax
2
at
(3.59)
reduces to
["
f"o(t)
.10.10
Gdxdt+c(x',t'),
(3.60)
having in mind the oxygen diffusion equation (1.57) and (3.57).
But
also,
integrating (3.59) by parts and using (3.56),
the
integral can
be
written
G--c-
dt-
[GcJ!~&>(")dx,
J.
t'
[
ac
aG]"="'<t)
11
o ax ax
,,=0
(3.61)
where t = to(x) is
the
inverse of x = xo(t),
the
path
of
the
moving bound-
ary in the
(x,
t) plane.
The
first integral in (3.61) and the contribution
from
t =
toCt)
in
the
second
one
vanish because
of
the two boundary
Integral-equation formulations
119
conditions (1.59)
and
(1.58).
On
inserting
the
initial condition (1.60) and
then equating (3.61)
and
(3.60) we find
1
,'
1"0(&)
1
i1
C(X',
t') = - G(x, x', t
'
- t)
dx
dt
+2:
G(x, x', t
'
)(I-
x f dx.
It
follows from (3.56) and (3.57)
that
1
~
a2G
G(x, x', t') =
[5(x-x
'
)5(t-
t
'
)+-2
(x, x',
t'-
t)]
dt,
ax
and
so (3.62) becomes
i
"
r
"0(')
C(X',
t
l
)=!(1_X
'
)2-
J
o
G(x, x',
t'_t)
dx
dt
(3.62)
(3.63)
(3.64)
Integrating
the
last integral of (3.64) twice by parts with respect
to
x and
using (3.58) and (3.56) we obtain
1
,'
C(X
I
,t
l
)=!(I-x
'
f-
G(O,x',t'-t)dt+
l
"i
1
+ G(x, x',
t'-
t)
dx
dt.
"0(')
(3.65)
We
can now substitute
the
chosen form of G(x, x',
t'_t)
from (3.55)
and
obtain
C(X',
tl)=!(I-XI)2-2(:Yexp(
-
:~~)+
+X'
erfC~~tl)+
R(x',
t'), (3.66)
where
1
r"{
[X
(t)-X
l
] [
I-x'
]
R(x
'
, t')=2:
J
o
erfc
2(tl-t)!
-erfc
2(t
'
-t)!
+
[
XO(t)+x
l
] [
l+x'
]}
+erfc
2(t'-
t)!
-erfc
2(t'-
t)!
dt.
(3.67)
The
first integral in (3.65) leads
to
the
second
and
third terms in (3.66),
after suitable integration by parts,
and
the
second integral in (3.65) leads
to
the
quantity
R(x
'
, t').
If
now we were
to
let
x'
in (3.66) approach
XO(t')
and remember from
120
Analytical solutions
the first of (1.59) that c(xo(t'), t') = 0 we would obtain an integral equation
for the function
x = xo(t). Hansen and Hougaard (1974) point out,
however, that alternative integral equations can
be
obtained by differen-
tiating (3.66) with respect to
x'
or
t and in each case putting
x'
= xo(t').
They decide
to
use the first of these forms as being simpler than
the
second
or
than (3.66). Writing xo(t) , xo(t') as x,
x'
respectively for
convenience they obtain
x'
=
1-2
erfC~~t')+
(1-
x')erf(~v:'')+
(1
+
x')erfc(~::'')
-2
(£y{
1-ex
p
[
(1
~t~')2]+exp[
-(14~/X')2]}
1
I"
1 {
[(x-
X')2]
[ (x
+X')2]}
V'" 0 (t' - t)! exp - 4(t' -
t)
-1-exp
4(t' _
t)
dt. (3.68)
Once the function
x = xo(t) has been determined from this integral
equation, c can
be
evaluated from (3.66) everywhere in the range
O<x<xo(t)
where c(x,
t»O.
Hansen and Hougaard (1974) adopt a
simple iterative numerical procedure
to
find x = xo(t) from (3.68). Writing
~
=
ih,
where h
is
a suitable time step, they evaluate the integral by a
standard numerical procedure, actually Simpson's formula coupled with
Newton's three-eighths rule. They express (3.68) as
h = N({;,
~)/D({;,
~)
(3.69)
where
h =
1-
xo(ih), and, writing generally
t;
=
1-
xo(jh), i
,;;;j,;;;
i,
i {
U)2
2}
N({;
~)=~
L
~
exp(-Vi-t;
)_1_exp(-(2-
h
-t;»)
"
V",
j=l
(~-~)'1
'
4(~
-~)
4(~
-~)
(
1-h)
(2-
h)
(~)!
+2erfc
2V~'
-2erfc
2Vt;'
+2
'"
x { 1 + exp
(-(24~
{;)2) - exp (
~~)
}.
(3.70)
D({;,
t;)=
1+erf~~t;)-erfc(~~:).
(3.71)
The
coefficients
CXj,
of
course, depend
on
the
quadrature formulae
adopted.
The
iterative procedure based on (3.69) is'
fl
n
+
1
)
=
N(fi,
t;)/D(fi, t;),
(3.72)
and the starting value
flO)
is
obtained by quadratic extrapolation from the
values
h-l>
h-2'
h-3
calculated at the three immediately preceding times.
Crank and Gupta (1972a) derived a formula for
c(x', t') for small times
Integral-equation formulations
121
TABLE
3.1
Position
of
the moving
boundary
as
a
function
of
time
Time x Time
x
0.000 1.00000 0.1800
0.50109
0.0200 1.00000 0.1900
0.34537
0.0400
0.99918 0.1950
0.20652
0.0600
0.99180 0.1955 0.18708
0.0800
0.97155 0.1960
0.16266
0.1000
0.93501 0.1965
0.13284
0.1200
0.87916
0.1970 0.09175
0.1400 0.79891 0.1972 0.06708
0.1600
0.68337
0.1974
a
0.00000
a Estimated
by
extrapolation.
which
is
identical with the solution (3.66) above except for
the
term
R
=
R(x',
t'). Hansen and Hougaard (1974) used the first three terms in
(3.66) to compute
XO(ti)
for
t;
<0.02.
For
larger times they used (3.72).
Finally, they computed
c(x, t) for various values of x and t from (3.66)
and showed that the results are in good agreement with those obtained
from the small-time solution of Crank and Gupta (1972a) given by (3.66)
without the term
R(x',
t'). Positions of the moving boundary are repro-
duced from Hansen and Hougaard's paper (1974) in Table 3.1, and Table
3.2 shows a selection of their values of
c(x, t).
These results are compared with those of other authors in Tables 4.1,
4.2, and 6.6-8.
TABLE
3.2
Values
of
10
6
c
x
0.0 0.1 0.2
0.3
0.4 0.5
0.6
0.7 0.8
0.9
0.000 500000
405000
320000
245000
180000 125000
80000
4500020000
5000
0.002
449537 401927
319973
245000
180000
125000
80000
4500020000
5000
0.010
387162
365072
309949
243275
179804
124986
79999
4500020000
5000
0.050 247687
240175
218841
186952 148990 109634
72961
42029
18856
4629
0.100
143177
139294
128082
110787
89295
65892
43018
23059
8232
603
0.150
63087
60845
54503
44503
32353
19583
8251
1005
0
0
0.160
48823
46840
41136
32434
21927
11304
2890
0
0
0
0.180
21781
20.287
16066
9942
3523
0 0
0
0 0
0.190
9021
7817
4578
799
0 0
0
0
0
0
0.195
2884
1914
32
0
0 0 0 0
0 0