Crank J. Free and Moving Boundary Problems
Подождите немного. Документ загружается.
132
Analytical solutions
For
an isothennal condition
on
x = t, so
that
(3.122) are replaced by
K
aU
1
=_H
1
ax
'
x=O;
uz=-V,
um=O,
(3.128)
Goodman and Shea (1960) derived the following set of ordinary differen-
tial equations
(3.129)
(3.130)
(J
=_
V(t_
)_
(t-s)2d(Jz
z 2 s 12kz
dt'
(3.131)
s(t~=O,
Vt
He-
(Jl(t~
=0,
(Jz(t~=
-'3-
6K
.
(3.132)
Poots (1962b) applied the heat-balance integral method to single-phase
melting problems in the semi-infinite region, the cylinder, and the sphere,
by using a two-parameter quadratic profile.
He
also developed a short-
time series solution. Lardner and Pohle (1961) also considered cylindrical
problems.
Goodman (1961, 1964)
and
Boley and Estenssoro (1977) suggested
ways of incorporating variable heat properties in heat-balance integral
methods, though not specifically in relation to problems of phase change.
Imber and Huang (1973) also studied temperature-dependent thennal
properties and when the dependence is linear they concluded
that
the
freezing of a semi-infinite liquid is well approximated by
the
constant
heat-properties solution when the properties
at
the fusion temperature
are used. Oearly, the mathematical manipulation in heat-balance
methods for anything
other
than relatively simple problems can
be
very
lengthy. The choice
of
a satisfactory approximation to the temperature
profile
is
acknowledged to
be
a major difficulty in the heat-balance
approach.
The
use of a higher-order polynomial, for example, does not
necessarily improve the accuracy of the solution. Goodman himself
(1%4), Bell (1978), Langford (1973), and others quote illustrations.
Boley (1973) suggested that improvement
is
obtained when approximate
solutions, such
as
the simple heat-balance method produces, are superim-
posed by using convolution
or
Duhamel integrals,
or
by the method of
images. Boley and Estenssoro (1977) likewise introduced the idea
of
'multiple penetration depths' and expressed the temperature as ,a series,
each
tenn
of which was taken to vanish beyond an unknown penetration
Integral-equation formulations
133
depth below the surface. Thus for a half-space, initially at temperature
zero and with its surface raised to a constant temperature
To,
the
temperature
T(x,
t) is written
as
where
T(x, t)/To = I d;x
i
1 -
~
H(x,
qi)
m-l
(
)2
i=O
qi
H(x,
X;)
=
{1,
x <X;,
0,
x>x;,
and d
o
=l,
d;=d;(t),
qi=qi(t)
for
i;;;.:O.
It
is assumed that
qi(t»O,i=
0, 1,
...
, m
-1.
Integration of the weighted-moments expression
1
00
(
iPT
aB
k---
x"dx=O
o
ax
2
at '
n
=0,1,2,
...
,
2(m-1);
using T =
T(x,
t) from above leads
to
(2m
-1)
equations
to
be
solved for
the
(2m
-1)
unknown
d;
and
qi.
Solutions were given for the cases m = 1,
2, 3 and found
to
provide improved approximations to
the
exact solution
as
more penetration depths were included. No numerical solutions were
obtained for variable thermal properties but the fundamental relation-
ships were
set
out.
3.5.5. Heat-balance method with spatial subdivision
Noble (1975) suggested that the most promising way of systematically
improving accuracy would
be
by repeated spatial subdivision, using quad-
ratic profiles in each sub-region, as in finite-element treatments. Bell
(1978) solved the single-phase melting problem by adopting a modified
form of Noble's suggestion in which equal subdivision of the dependent
variable, temperature, is used rather than the more usual equal division of
the space variable.
The
solution
of
a system of differential equations
provides the positions of each isotherm identified by the subdivision,
including automatically
that
of the melting boundary, at successive times.
Bell (1978) used the equations
au=
k a2u
at
ax
2
'
O<x<so(t),
t>O,
(3.133)
u(x,
t) =
Uo,
so(t)~x,
t;;;.:O,
(3.134)
u(O, t) =
u.,
t>O,
(3.135)
_Kau=Lp
ds
ax
dt'
x = so(t),
t>O.
(3.136)
The
temperature range
Uo
to
U.
is
subdivided into n equal intervals as in
134
Analytical solutions
U
Us
S;
So
x
FIG. 3.1. SUb-division of region
Fig. 3.1.
The
depth of penetration
of
the isotherm
u;
at time t
is
denoted
by
Si(t)
where
i
=0,1,
...
, n,
and
Un
=
u..
Integration
of
(3.133) over the sub-region
S;+1
(t):S;;;
x:S;;;
s;
(t)
gives
~
[is.
u(e,
t)
de]
-
u;
ds
i
+
U;+1
dS
i
+1
= k
(au)
_ k
(au)
,
dt
"+1
dt dt
ax i ax
i+1
more conveniently written as
d
d
[IS,
U(e,t)d
e
+U;+lS;+l-U;
S
i]=-k(au).
+k(au).'
(3.137)
t
"+1
ax
1+1
ax ,
which, remembering (3.136), takes the special form
:t
[fa
u(e,
t)
de
+
U1
S
1-(UO-
LIC)So]
=
-k(!;\,
Sl(t):s;;;X
:s;;;so(t).
(3.137a)
Clearly, (3.137) expresses the heat balance of the sub-region
S;+1(t):s;;;x:s;;;
Si(t).
Considering each sub-region
(S;+1,
Si),
i = 0, 1,
...
, n
-1,
in turn, n
equations
of
the form (3.137) are produced in terms of the n unknowns Si
(i
= 0, 1,
...
, n
-1)
and the corresponding known temperatures
U;,
recal-
ling that
s" = 0.
The
heat-balance equation commonly used
in
the Good-
man integral method
is
simply the sum of these n equations.
The
temperature u(x, t)
is
. approximated by a series of profiles VI which'
Integral-equation formulations
135
simultaneously satisfy the n equations (3.137), i =
0,1,
...
, n
-1.
Thus
Vi(X,
t)
=
u(x,
t) in the sub-region
S;+I:5;;X
:5;;S;,
i =
0,1,
...
, n
-1,
and
hence
Vi
(Si)
=
U;
and
Vi
(S;+l)
=
U;+I'
If
Vi
is
chosen to
be
a quadratic
function in
x then one extra condition
is
needed
to
specify
Vi
in each
sub-region. These can
be
the conditions that express the continuity of
heat flux between each subdivision and its neighbours, together with the
special condition (3.136)
at
x = So for the end subdivision.
For
conveni-
ence, let
so that a suitable quadratic fonn for
Vi
is
Vi
(Si+l
+
Yi)
=
U;+1
-
Yi(U;+1
-
U;)/.f;
+
a;Yi(1-
yJ(i),
(3.138)
where the coefficient a;
is
detennined from the condition for continuity of
heat flux
(aVJax)
..
=
(a1J;_I/aX)
..
,
The
recurrence relation
u.-uo
(1
1)
a;
==--
---
-a;-l>
n
(i-l
t;
(3.139)
(3.140)
follows immediately from (3.138) and (3.139), so that, once
ao
is
deter-
mined from (3.136) the a; can
be
generated.
Mer
differentiating (3.138)
and
putting i = 0,
Yo
= (0 to obtain (au/ax)
..
and
(a
2
u/ax
2
)
..
, use of the
condition (3.136), rewritten in the form
of
(3.110), namely
(3.141)
leads to the result
(3.142)
where
y = Ln/{(u. -
uo)c},
(3.143)
and the negative square root is chosen
to
fit
the physical fact that as the
latent heat
L becomes larger, the melting boundary moves more slowly
and the temperature distribution approaches a linear steady state. This
requires, from (3.138) that a;
~
0, i = 0, 1,
...
, n
-1,
and in particular
ao~O
as
L~oo
which
is
satisfied by (3.142)
and
(3.143).
Combination of (3.138), (3.142), and (3.137a) gives the equation
d {
n-l
}
6k
-d
(3-{3+6Y)(0+6(1+y)
L
~
='0
(l+{3),
t
/=1
1:'0
(3.144)
136
Analytical solutions
where
n-l
L
~=SI'
j=1
Similarly, after generating
the
as
from (3.140)
and
substituting (3.138)
into (3.137),
the
remaining
(n
-1)
equations
are
found
to
be:
~
[
2{
i-I
(_1)i-i-
1
(-
W-
1
(1
+
(3)}
n-l]
dt
(i
2
j~1
~
+
(0
+2(i
+6
i=~1
~
_
12k
{i-l
(_1)i-i-
1
(-
W-
1
(1
+
m}
-g--12k
2 L
e,
+
e.
,1~i~n-1,
(3.145)
~i
j=1
i 0
which together with (3.144) form a system
of
n first-order non-linear
differential equations for
the
n penetration parameters
(i
(t), 0
~
i ~
n-1.
Bell (1978) suggested
that
the
initial singularity in
the
equations can
be
dealt with
by
use
of
a series solution, e.g.
that
of
Poots
(1962b),
but
refers
his own numerical calculations
to
the
situation in which
the
penetration
of
the
melting interface is proportional
to
(kt)!. Solutions
of
the
form
(i
='\.;(kt)!,
A;
>0,
are
defined
by
reduced forms
of
(3.144) and (3.145),
namely
(3 -
(3
+
6y
)A~
+ 6(1 + Y
)Ao
{~:
"-i
}-12(1
+ m = 0 (3. 144a)
and
(
i-l
(_1)H-l
(-1Y-l(1
+
f.l))
A~
L +
,.,
+Af
•
j=1
"-i
2Ao
•
(
n-l
i-I
(-W-i-
1
4(-W-1(1+m)
+3A;
i=tl
"-i+8i~
"-i
+
Ao
-12=0,
1
~
i
~
n-1.
(3.145a)
Bell (1978) started a Newton-type iterative process
of
solution using
initial estimates
of
AI'
...
,An-I
taken
from
Goodman's
integral solution
(1964), which corresponds
to
the
case n = 1 in Bell's method. Table 3.3
compares Bell's results for
the
position
of
the
melting interface so/(kt)!,
for different subdivisions, with those for
the
exact solution (Carslaw
and
Jaeger
1959) for
the
case L/(u. -
uo)c
= 1.
TABLE
3.3
Position
of
melting interface/(kt)! i.e. so(t)/(kt)!; L/(u. -
uo)c
= 1
Exact
n
1
2
4
8
32
solution
1.32003
1.25452
1.24291 1 240 70 1.240 15 1.240 13
Integral-equation formulations
137
TABLE
3.4
Percentage
error
in position
of
melting
front/(kt)!,
i.e. in so(t)/(kt)!, compared with
exact solution. Two sub-regions n
= 2: differ-
ent
latent heat parameters; L/(u. - uo)c
L!(u,.-uo")c
% error in
so(t)!(kt)f
1
10
50
+ 1.16 +0.24 +0.05
100
+0.02
Table
3.4
shows results obtained when only two sub-regions
are
used
for different values
of
the
latent
heat
parameter L/(u. -
Uo)c.
This
method, based
on
equal subdivision of
the
dependent variable tempera-
ture, seems a promising
one
to
apply
to
other heat-flow problems with
or
without phase changes.
In
a
later
paper, Bell (1979) applied his spatial subdivision method
to
the
solidification
of
a liquid around a cylindrical pipe
but
did
not
assume a
similarity solution.
The
more
general equations were solved by a
Runge-
Kutta
method,
and
also by a simpler numerical integration algorithm,
starting from a short-time series based
on
Poots's (1962b) solution.
3.5.6. Multi-dimensional problems
Poots (1962a) extended
the
heat-balance integral method
to
study the
movement
of
a two-dimensional solidification front
in
a liquid, initially at
the
fusion temperature, contained
in
a uniform prism of square
c1:"oss-
section
of
side
a.
The
surface is maintained
at
a constant temperature
below
the
fusion
temperature
of
the
liquid. Constant thermal properties
are
assumed.
The
prism is
bounded
by
the
isothermal surface
at
temperature
To,
where
To
< T
F
,
the
fusion temperature
of
the
liquid.
At
any time, while
the
process
of
inward solidification is occurring, the
two-dimensional solidification front is defined
to
be
the
surface
Cp(x, y, t) = 0, an isothermal
on
which T = T
F
•
After
introducing
the
dimensionless variables
f) =
T-To
_
pLk
T
F
-
To
'
{3
-
K(T
F
-
To)
,
(3.146)
.f=s/a,
X=x/a,
Y=y/a,
138
Analytical solutions
the
definitive equations may
be
written
and
where
a
2
8 a
2
8 a8
aXZ+
aYZ=
aT'
8=0
on
Co(.~
Y)=
(XZ-1)(YZ-1)
=0,
8=1,
Cp(X,
Y,T)=O,
f
a8 dv
-dt=-/3-
Cp=O
all
dT'
Co=o
v(Cp,C
o
)=
II
dXdY.
Cp=O
The
initial condition is
8=1,
(3.147)
T>O}
(3.148)
(3.149)
(3.150)
In
these equations, t is
the
non-dimensional distance along Cp(x,
y,
t) = 0
in an anti-clockwise direction
and
II
is
the outward normal
to
this
solidification front which
is
assumed
to
have
the
shape
Cp(x,
Y,T)=(XZ-1)(YZ-1)-e(T)=0
- (3.151)
where
e
is
an unknown function
of
time
T.
The
initial condition (3.150)
now becomes
e=O,
1'=0,
(3.152)
and (3.151) implies that
e = 1
at
the
instant solidification
is
completed.
The
heat balance integral is now obtained by integrating
both
sides
of
(3.147) over the solid phase bounded by
the
contours Co=O
and
Cp=O.
Application
of
Green's
theorem
to
the
result and use of (3.149) leads
to
the
relationship
(3.153)
A second integral
Co=o
/3
::
= I I
{~;Y+(:;Y
+8
::}dXdY
(3.154)
Cp=O
Approximate solutions
139
follows by applying Green's theorem
to
the divergence theorem
div(6 V6) = V6 .
V6+
6
~6
= V6 . V6 + 6
a6,
aT
on substituting from (3.147), and using (3.148) and (3.149). Poots
(1962a) completed the analysis for an assumed one-parameter tempera-
ture distribution
6 =
(XZ-1)(YZ-1)/e
(3.155)
and also for a two-parameter form
(XZ-1)(YZ-1)
(XZ-1)(YZ-1))2
6=
(l-g)+g
.
e e
(3.156)
Ivantsov (1947) showed that for an isothermal phase-change surface
which moved with a constant velocity, the shape of a paraboloid of
revolution satisfied
the
differential equation and boundary conditions.
Horvay and Cabn (1961) extended
the
same approach
to
other assumed
surface shapes.
Riley and Duck (1971) considered the three-dimensional freezing of a
cuboid.
3.6. Approximate solutions
3.6.1. Steady-state approximation and improvements
In
Neumann's similarity solution discussed in §3.2, the equation (3.14)
describing the motion of the solidification boundary in the single-phase
problem was seen
to
reduce
to
(3.15) for A small, which corresponds to a
latent heat L that
is
large compared with the specific heat capacity of the
ice.
Stefan
(1889a) and sundry authors later obtained the result (3.15) by
using the physical approximation that for large L, implying a slowly
moving solidification boundary, the temperature distribution at any in-
stant behind
the
moving boundary
is
the steady-state distribution which
would
be
set up
if
the boundary were
to
be fixed in its position at that
instant.
In
effect, the heat
flow
is
rapid enough to maintain this pseudo
'steady state'.
Thus, putting
aUl/at
= 0 in (3.1), the steady-state solution which satisfies
(3.4)
and
the first of (3.6)
is
(3.157)
140
Analytical solutions
Substitution of (3.157) in the moving-boundary condition, -
Kl
au1/ax =
Lp
ds/dt, for the single-phase problem yields
S2
=
2Kl
U1t/Lp (3.158)
which
is
equivalent to (3.15) remembering that in (3.14) A = a/(2kI) =
!S/(klt)~
from (3.10), where
kl
==
K1/C1P.
Crank (1975, p. 310) used
the
steady-state approximation
to
solve the
diffusion version of the one-phase melting problem for a surface bound-
ary condition
e au/at = K au/ax,
x=o,
t>o,
(3.159)
with e and K constants. Crank (1957b) gave the corresponding approxi-
mations for the cylinder and the sphere together with estimates of
accuracy based
on
numerical calculations.
The essential simplifying feature of the above method
is
that
the
assumption of the steady-state temperature profile avoids the need to
solve the heat-flow equation. Various authors have exploited this advan-
tage by assuming other more sophisticated profiles, of which the steady-
state solution
is
usually the leading
term
in a series. The higher terms may
therefore provide an estimate of the errors associated with the steady-
state assumption.
Thus, Solomon
(1978) described an approximate solution, following
Megerlin
(1968), for a constant surface heat
flux.
The method follows the
same pattern as that of Goodman
(§3.5.4) but no attempt
is
made
to
preserve the heat balance. Take
the
problem defined by the following:
Consider
o<x<s(t),
au/ax
=-1,
x
=0,
au/ax
==
-ds/dt,
u
=0,
u=o,
x
>s(t);
t>o,
x = s(t),
s(O)
=0.
t>o,
(3.160)
(3.161)
(3.162)
(3.163)
(3.164)
where
at>
az,
s are functions of t to be determined and the second of
(3.162)
is
automatically satisfied. Using (3.162) and its alternative form
(see equation
(3.110)) coupled with differentiation of (3.164) it im-
mediately follows that
ds 1 (dS)2 2
u=
--(x-s)+-
-
(x-s).
dt
2
dt
(3.165)
Approximate solutions
141
Substitution of (3.161), after differentiation of (3.165), leads to
ds
-1
+.J(l
+4s)
-=
dt
2s
the negative root being discarded on physical grounds, and finally to the
solution
(3.166)
We note that, in this problem, the steady-state approximation
is
u = s - x
and
dsldt=
1, which corresponds to the first term in (3.165). For this
term, therefore, a heat balance over the whole solid region is preserved.
Solomon (1978) compared his results computed from (3.166) with
those obtained by Rose (1960) for the same problem based on a weak
solution. Table 3.5
is
an extract from Solomon's table (1978).
Solomon, Wilson, and Alexiades (1981) showed that the solution to the
Stefan problem with a convective boundary condition tends to the quasi-
steady approximation as the specific heat tends to zero.
Pekeris and Slichter (1939) study the formation of ice around a cold
cylinder embedded in water initially
at
O°C,
starting from an assumed
temperature profile given by
u
=A+B
In
x+
(~)H(A
-B)x
2
+;lBx
2
lnx}
+(::){1~:
(2.4-38)+!
x4Inx}+
...
, (3.167)
where A,
B are functions of time, B = dB/dt etc., a is the radius of the
cylinder, and
x =
rIa.
They finally arrive
at
the solution, for y = {s(t)la}2,
TABLE
3.5
Comparison
of
t(s) from (3.166)
and Rose (1960)
set)
t[Eqn (3.166)]
t(Rose)a
0 0
0
0.4 0.450
0.466
0.8 1.030
1.034
1.2 1.690
1.681
1.6 2.440
2.394
a The agreement
is
stated to
be
within
the limits of accuracy of Rose's values.
(3.168)