Crank J. Free and Moving Boundary Problems
Подождите немного. Документ загружается.
212
Front-fixing methods
where u
denotes
temperature
and
constant thermal properties
are
as-
sumed.
The
radial coordinate,
r,
is measured from
an
origin
that
remains
fixed as time,
t,
changes. By analogy with (5.59)
an
IMM form of (5.111)
is
:;
=
~::)
/
(:~r
-~-
:2
(::)
::~.
(5.112)
In
general, flow
of
heat
in two space dimensions
can
be
represented
by
a family
of
isotherms
and
an
associated family
of
flow lines, orthogonal
to
the
isotherms.
In
an
isotropic medium, any
point
on
an
isotherm moves
along
the
flow line normal
to
the
isotherm
at
that
point.
Heat
flow is
everywhere
normal
to
the
isotherms
and
never across flow lines.
By
confining attention
to
a small segment
of
an
isotherm for a
short
interval
of
time,
we
can
regard
the isotherm element as
part
of a cylindrical
system
and
identify
the
coordinate
r,
in (5.112), as
the
local radius
of
curvature
of
the
isotherm measured from
the
local
centre
of
curvature
assumed fixed in its position
at
time
t.
Equation
(5.112) yields
the
velocity
iJr/iJt
of
the
selected
element
along
the
normal
to
itself.
Because
the
general system is distorting
and
rotating,
both
the
centre
of
curvature for
the
element
and
the
curvature itself
may
change with time
as well as from
point
to
point
in
the
system.
The
flow lines will
not
strictly
be
radial lines
of
constant
6,
and
the
local isotherms will
not
be
exactly
concentric circular arcs. But, because r in
(5.112) is chosen
to
be
along
the
local normal,
the
term
iJ
2
u/iJ6
2
will
be
small in general,
but
non-zero.
Crank
and
Crowley (1978) approximated
it
to
zero
and
calculated
the
movement
of
a
point
on
an
isotherm along its normal, in a small time
at,
by
solving
the
equation
:;
=
(:::)
/
~:r
-~
.
(5.113)
They describe
an
explicit finite-difference
method
based
on
three
typical
isotherms (Fig. 5.6)
on
which
the
temperatures
are
(j-1)
5u,
j5u,
and
FIG.
5.6. Sketch showing relative positions of isotherms
Isotherm migration method
213
(j
+ 1) 5u and
A,
B, C are three points whose coordinates are known on
isotherm
j 5u.
The
points G and F are the intersections of the radius
rm
with the chords approximating
the
isotherms (j
-1)
5u
and
(j
+ 1)
5u.
Let
~
=
Irml
be
the
distance of B from the centre of curvature
of
the arc
ABC
and
let
n!, and
n;;'
denote the distances of F and G from this same
centre respectively. Then the usual approximations are
an)
= n!, -
n;;'
a
2
n)
= n!, -
2~
+
n;;'
(5.114)
au
B
25u'
au
2
B (5u? .
In
the
interests of accuracy, it
is
preferable to calculate the differences
n!, -
n;;',
n;:'
-
~
and
~
-
n;;'
directly from the coordinates of the points
B, F,
and
G which are known at time
t.
Thus, denoting by
aHj,m
the
movement of
the
point m
on
the
isotherm j 5u along the normal in the
small time interval from t to t
+
at,
Crank and Crowley (1978) replaced
(5.113) explicitly by
at
~
(5.115)
If
the coordinates of the point m
on
isotherm j
5u
at time t = i
at
are
denoted by
xtm,
ytm,
then the new coordinates
at
time
(i
+
1)
at
are
H1
_ ; +
A;
O;}
Xj,m
-
Xj,m
anj,m
cos
j,m,
y;'~1
=
ytm
+
antm
sin
O;'m,
(5.116)
where
8;'m
denotes
the
direction of the normal
to
the isotherm at B
measured from the
x-axis,
as
in Fig. 5.7.
It
remains
to
calculate
O;'m
from
the coordinates
of
the
points A, B, C, and also
the
length of the radius of
curvature,
~,
at B (x".,
Ym).
Figure 5.7 shows a section of
an
isotherm approximated by a circular
arc
ABC
which
is
concave downwards and to which
is
assigned a positive
FIG. 5.7. Sketch illustrating the geometrical procedure used
to
determine the normal
to
each isotherm
214
Front-fixing methods
curvature.
The
tangents
at
the
mid-points
P,
Q
of
each
of
the arcs
AB,
BC
are
parallel
to
the
corresponding chords
AB,
BC.
Thus
the
change in
the
direction
of
the
tangent along arc
PBQ
is given
by
I/Im
-I/Im+t.
where
the
t/ls
are
the
angles
made
by
the
perpendicular bisectors
of
the
chords
AB
and
BC
respectively with
the
x-axis
in
Fig. 5.7.
The
arc length
PBQ
may
be
approximated by !(s",+1 + s",),
where
s",
denotes
the
length
of
the
chord
AB,
labelled chord m.
Then,
the
radius
of
curvature
rm
at
B
(x,..,
Ym)
may
be
written
~
S)
1
(s",
+
s",+1)
~-
-
--
-
t/I
m - 2
t/lm
-
t/lm+l
.
(5.117)
By regarding
BQ
as a circular arc
of
radius
Irml
=
~
and
length
!s",+1>
the
angle
marked
'Ym
on
Fig. 5.7
may
be
approximated as
is",+l/~.
Hence
the
angle 8:'
m
in (5.116), illustrated by 8
m
in Fig. 5.7, can
be
calculated
from
(5.118)
The
angles
t/I
and
distances s
are
readily
computed
from
the
coordinates
of
the
points
A,
B,
C.
If
the
initial
data
in
a two-dimensional heat-flow problem
are
given as
the
coordinates
of
a
set
of
points
on
each
of
a
number
of
isotherms,
the
relationships (5.115-118) provide
the
basis
of
a numerical algorithm
to
advance
each
point
in
a
succ~ion
of
small time steps,
At.
In
a Stefan
problem
in which
the
phase-change boundary is itself
an
isothermal surface, its motion is readily calculated
by
the
IMM
form
of
the
usual Stefan condition, e.g.
~;=-~~:rl,
(5.119)
in
a one-phase problem, which, suitably discretized, is used in place
of
(5.115).
Crank
and
Crowley (1978) solved
the
problem
of
the
solidification
of
a
square prism
of
fluid
set
out
in §5.4.2(i).
They
discussed questions
of
accuracy
and
stability heuristically
and
obtained
numerical results which
are
in good
agreement
with those obtained
by
Crowley (1978) (see §§6.2,
6.2.4)
and
by
Elliott
(1976) as shown in Tables 6.12, 13.
The
initial
isotherm positions
needed
to
start
the
IMM
calculation were
obtained
by
the
enthalpy method.
Because
of
symmetry
about
the
diagonal Y = x it is sufficient
to
work
in
the
triangular region 0 ~ x ~
Y,
0
~
Y ~ 1.
The
end
points
on
the
bound-
aries x = 0
and
Y = x
need
special consideration. A circle, with its
centre
on
the
axis
at
(0,
Yo)
or
on
the
diagonal
at
(a,
a)
as appropriate, is fitted
through
the
end
point
and
the
next
point
inside
the
region
on
each
Isotherm migration method
215
isotherm, in
order
to
find
the
curvatures. The angles 8
i
.
O
=
'Tr/2
on
the
y-axis
and
8
i
.M,
say, =
'Tr/4
on
the diagonal are already known. Thus
on
the
axis, x
2
+ (y -
Yo)2
=
r2,
and we find froro (0,
Yi,O),
(Xj,1>
Yi,l)
that
r = [(Yi,o-
Yi,l)2+
xT.l]/2(Yi,o-
Yi,l)'
On
the
diagonal
the
use of
(x
-
af+
(y-a)2=
r2
yields
r
=
{(Xj.M
-
Xj.M_l)2
+
(Yi.M
-
Yi,M_l)2}/2(Xj.M
+
Yi,M
-
Xj,M-l
-
Yi.M-l),
and
one
or
other
of these expressions for r replaces (S.117).
Crank
and
Crowley (1979) described a corresponding implicit method
for tracking isotherms along orthogonal flow lines, using a linearized form
of equation (S.113).
In
general,
the
implicit scheme offered
the
usual
advantage over
the
explicit method in enabling the use of longer time
steps. However, they choose
to
illustrate their implicit scheme by consid-
ering
the
problem of Bonnerot
and
Jamet (1977) specified in equations
(4.24-8). This problem is a more exacting test of Crank
and
Crowley's
method for two important reasons. Firstly, the radius of curvature is
positive
at
the
left-hand
end
of each isotherm, and negative as the
isotherm approaches
x = 1 (Fig. 4.4(a»). Therefore, proceeding in
the
direction of x increasing along any isotherm the curvature increases from
finite positive values, through infinity,
and
then from negative infinity
to
finite negative values,
the
central
part
of each isotherm being virtually
straight.
In such a situation, care is necessary
to
ensure that
the
finite-
difference equations are written correctly when
rl,m
<0
and
"'.t,m
= -ri,m'
The
second difficulty is associated with
the
initial temperature distribu-
tion.
On
the
left boundary, x = 0,
iJr/iJt
in (S.113) can
be
identified with
iJy/iJt
since
on
x = 0
the
isotherms are concave downwards
and
the
centres
of curvature are
on
x =
O.
When
the
necessary derivatives are obtained
from
the
second of
the
conditions (4.27) and substituted into equation
(S.113) we find that
on
x=O,
iJy/iJt=-'Tr2(1-u)~0,
O<u<1.
However,
from
the
first of (4.28)
on
the
phase-change boundary, u
=0,
we have
iJy/iJt=1!3A
on
x=O.
Thus,
at
t=O, all the isotherms except
u=O
initially move in
the
direction of Y decreasing
at
the
left boundary x = 0,
whereas
the
phase-change isotherm, u = 0, moves upwards. In fact, there
is a discontinuity in
iJy/iJt
at
t = 0
on
x = 0, and it is
not
surprising that
the
finite-difference solution of
Crank
and Crowley (1979) showed a loss of
accuracy there. Accuracy can
be
improved by adding an extra isotherm
with the temperature
~
Bu.
On
the right boundary, x = 1, all isotherms move upwards initially, since
iJy/iJt
=
'Tr
2
(1-
u), 0 <
u.,;;;
1, x = 1,
and
for u = 0,
iJy/iJt
=
1/A,
x =
1.
These
statements are physically consistent with
the
initial temperature distribu-
tion in (4.27).
The
negative temperature gradient from left
to
right along
any line of constant y produces a corresponding sideways heat
flow
which
causes
the
isotherms to move downward
on
x = 0 and upward
on
x =
1.
216 Front-fixing methods
Crank and Crowley (1979) presented numerical results. Some compari-
sons with those of other workers for the positions of the ends of
the
interface are reproduced in Table 5.2.
TABLE
5.2
Comparison
of
positions
of
ends
of
interface
for Bonnerot and Jamet's
(1977) problem
specified in equations (4.24-28)
IMM
JBa
RMP Enthalpy" implicit
d
On
x=O
0
3
3 3
3
1 3.021 3.015 3.031 3.073
2
3.068 3.054 3.051 3.161
On
x=1
0
1
1 1
1
1
2.118 2.095 2.124 2.132
2
2.610
2.585
2.566 2.591
a
JB
= Bonnnerot and Jamet (1977).
b
RMF
= Furzeland, R. M. (1977b).
IMM
explicit
e
3
3.057
3.131
1
2.133
2.594
C Enthalpy results obtained by method
of
Crowley
(1978).
dIMM
implicit=Oank
and Crowley (1979).
e IMM explicit results obtained by method of Crank
and Crowley (1978).
6.
Fixed-domain methods
6.1. Introduction
THE
previous three chapters have described what are essentially 'front-
tracking' methods. Even though, in Chapter
5,
the
phase-change bound-
ary was 'fixed'
by
a change of variable which simplified
the
numerical
work considerably, it was still necessary
to
satisfy
the
Stefan
or
similar
derivative condition
on
that boundary. Furthermore, it may sometimes
be
difficult
or
even impossible
to
track
the
moving boundary directly
if
it
does
not
move smoothly
or
monotonically with time. Quite possibly, it
may not do so in more than
one
space dimension.
The
moving boundary
may have sharp peaks,
or
double back,
or
it may even disappear.
The
possibility, therefore,
of
reformulating the problem in such a way that
the
Stefan condition is implicitly
bound
up in a new form of
the
equations,
which applies over
the
whole of a fixed domain, is an attractive one.
The
position of
the
moving boundary appears, a posteriori, as
one
feature of
the
solution.
One
way of accomplishing such a reformulation is
to
introduce an enthalpy
or
total heat function.
6.2. Enthalpy
method
The
use of enthalpy was proposed by Eyres et
al.
(1946) and
later
by
Price
and
Slack (1954) and Albasiny (1956). Rose (1960) introduced
the
concept of a weak solution used by Lax (1954)
to
obtain solutions of
hyperbolic equations involving shocks. Oleinik (1960) and
Kamenomostskaja (1961) provided theoretical justification.
The
proce-
dure
is
to
introduce an enthalpy function,
H(u),
which is
the
total heat
content, i.e.
the
sum of
the
specific
or
sensible heat and
the
latent heat
required for a phase change.
The
heat jump
at
the
phase-change bound-
ary is incorporated in
the
definition of
H(u),
which is
H(u)
=
1:
p(6)c(6)d6,
u<u
m
,
H(u)
=
1:
p(6)c(6)d6+pL,
u>u
m
,
(6.1)
1:
p(6)c(6)
d6~H(u)~
1:
p(6)c(6)
d6+pL,
218
Fixed-domain methods
H
o u
FIG. 6.1. Enthalpy function
In
(6.1)
Uo
is some fixed temperature less than
u.n,
the
melting tempera-
ture: density
p
and
specific heat c may both
be
functions of temperature
u.
Budak et
al.
(1965), Meyer (1973), Bonacina et
al.
(1973) incorporated
the
three expressions in (6.1) into the single relationship
H(u)
=
lU
{p(6)c(6)+Lp(6)8(6-unJ}d6,
Uo
(6.2)
where 8 is
the
Dirac impulse function.
The
form
of
H(
u) is shown
graphically in Fig. 6.1.
Consider
the
two-phase, multi-dimensional problem
i=1,2,
[K
au/anN =
-pLv,.,
x = s(t),
(6.3)
(6.4)
where
i = 1, 2
denote
the two phases separated by
the
moving phase-
change boundary
x =
s(t):
U;(x, t), i = 1, 2, denote
the
temperatures u in
the
two phases respectively: and
1:;,
Pi,
K;.
may all
be
functions of
U;,
x,
t,
where x
is
in general a three-dimensional space vector;
v,.
is the velocity
of
the
moving boundary in the direction of
the
normal,
n.
For
variable
properties,
c(u), p(u),
it
follows from
the
first two expressions of (6.1)
that
aH/at =
p(u)c(u)
au/at,
(6.5)
and
hence the problem of (6.3) and (6.4) can
be
:reformulated over
the
whole fixed domain occupied by the two phases together, except where
u
::::
Urn, as
the
single equation
aH(u)/at
=
V.
(KVu),
(6.6)
with associated initial and fixed boundary conditions.
Enthalpy method
219
Equation
(6.6) is not obviously meaningful at u =
Urn
since
the
enthalpy
H(u)
has a jump discontinuity
at
u =
Urn
of
magnitude [H(unJJi =
Lp
and
dH/du
is infinite there. Shamsundar and Sparrow (1975), by considering
the
energy balance of a control volume V with surface
area
A,
set up an
integral form of
the
enthalpy equation, which includes fluid motion due to
density changes
or
convection, i.e .
.:!.
r H dV+J
Hv·
dA=
J
KVu·
dA,
dd
v
A A
valid throughout
the
domain.
The
second term
on
the
left-hand side
accounts for any fluid motion. By applying this integral form
to
an
element of volume which does
not
include the phase-change boundary
and
then
to
one
through which
the
interface passes, Shamsundar and
Sparrow used
the
divergence theorem
to
confirm that
the
integral form is
equivalent
to
(6.3) and (6.4) and
to
(6.6) when
v=O.
In a one-phase problem, it follows that if the liquid is uniformly
at
the
melting temperature initially, we have (Shamsundar
and
Sparrow 1976)
:ti
H\dV+ L
H\v·dA=O,
where
H\
is
the
enthalpy of
the
liquid, and so
on
subtraction from
the
integral enthalpy equation above
.:!.i(H-HJdV+J
(H-H\)V.dA=J
KVu·dA.
dt
A A
In
the
solid phase v = 0 and so
an
enthalpy equation which takes account
of density changes in a one-phase problem is
:tL(H-HJdV=
L
KVu·dA,
which avoids explicit use of
the
fluid velocity.
For
constant thermal properties
Cl
and
C2 over each phase, constant p
over both phases,
and
urn=O,
the
enthalpy function defined in (6.1)
becomes
H(u)
= PC1U,
u<O;
H(u)
= PC2U + pL, u
>0;
O:;o;H:;O;pL,
u =
O.
(6.7)
If
K =
K(u)
in (6.6)
the
equation can
be
made linear by what
is
some-
times called
the
Kirchhoff transformation (Eyres et al., 1946; Albasiny,
1956):
(6.8)
220
Fixed-domain methods
which reduces (6.6)
to
(6.9)
This transfonnation can sometimes simplify
the
nU!flerical work in
other
methods too.
For
constant conductivities,
Kl
and
K2
in
the
two phases, (6.8)
be-
comes v = K;u, i = 1, 2,
and
H(
v) is given by
H(V)=PC1V/Kl>
v<O;
H(v)
= PC2v/K2+pL,
v>O;
O,,;;;;H(v)";;;;pL,
v
=0.
(6.10)
One-phase problems can
be
treated
by
extending
the
definition
of
u
to
a fixed, two-phase domain in which
the
usual two-phase equations
hold
but
c(u)=p(u)=O
in
the
second
phase
(Budak
et
al.
1965). Effectively,
(6.6) is then
to
be
solved only in
the
first phase.
6.2.1. Weak solutions
We
have
noted
that
H(u),
as defined in (6.1), is a discontinuous
function
of
u
at
u = Urn
and
this presents difficulties in interpreting
equation (6.6)
and
its classical solutions
at
u = Urn.
In
such situations,
the
usual practice is
to
seek
a 'weak' solution
that
satisfies a suitably integ-
rated
fonn
of
the
equation in which derivatives
of
Hand
u
do
not
appear
(e.g. Williams, 1980,
pp.
51, 142).
Weak
solutions for Stefan problems
were introduced
by
Rose
(1960), Oleinik (1960),
and
Kamenomostskaja
(1961),
and
have
been
extended
and
analysed
by
many
authors including
Friedman (1968), Brezis (1970),
and
Lions (1969).
Many
references
are
given in
the
survey
by
Niezg6dk:a (1983).
The
derivation
of
a
weak
solution can
be
illustrated
by
considering a
one-dimensional problem.
Equation
(6.6) becomes
aH=~{K(U)au},
0<x<1,
(6.11)
at
ax ax
and
we
take
boundary
and
initial conditions (Fox, 1979)
u(O, t) = gl(t),
u(1,
t) = g2(t),
u(x,
0) = uo(x),
(6.12)
and
consider
the
time
range
0,,;;;;
t,,;;;;
T.
We
first confine attention
to
this
problem in
the
absence
of
any phase change
so
that,
instead
of
the
definitions (6.1),
we
assume
the
higher-order derivatives
of
H and u
to
exist.
We
now multiply (6.11)
by
an
arbitrary
test
function q,(x, t) which
has
continuous first
and
second derivatives with respect
to
x
and
t
and
which vanishes
on
x=O
and
x=1.
Thus
from (6.11)
we
obtain
fllT[q,~
(Kau)_q,
a!fl
dx
dt=O,
J
o
ax ax
at
J
Enthalpy
method
221
and,
on
integrating by parts, the left-hand side becomes
I
T[
au]l
I11T
auact>
ct>K-
dt-
K--dxdt
ax
0
ax ax
-
f[ct>H]~
dx+
flTH::
dx
dt
I
T[
act>]
1
I11T
a (
act»
= -
uK
ax
0
dt
+ u
ax
K
ax
dx
dt
+
fH(uo)ct>(x,
0)
dx+
{11TH:: dx
dt,
and
so finally we obtain the required weak form
f
l1T{
a
(act»
act>}
u-
K-
+H(u)- dx
dt
o
ax
ax
at
I
T
act>
IT
act>
= K(g2)g2 - (1,
t)
dt-
K(gl)gl
- (0,
t)
dt
ax ax
-
fH(uo)ct>(X,
0) dx. (6.13)
We
now define a weak solution of a heat conduction problem with a
phase change, subject
to
the boundary and initial conditions (6.12),
to
be
a
pair
of bounded functions u and H, related by (6.1), for which the
integral form (6.13) is satisfied for all test functions
ct>
as specified.
Atthey
(1974) showed
that
the
weak formulation (6.13) includes the
usual Stefan
jump
condition (6.4) across a phase-change surface. Thus,
we assume that, as in the classical Stefan formulation, there are two
regions
Rl
and R2 separated by a curve x = s(t) in Fig. 6.2, and in each
region the one-dimensional form
of
the
heat conduction equation,
pc
au/at
=a(K
au/ax)/ax,
holds.
On
integrating by parts over
Rl
as in the
f=Tr--~-------'---~
x-O
x-I
FIG. 6.2.
x-t
plane for weak solution containing jump discontinuities
on
x =
set)