Crank J. Free and Moving Boundary Problems

222

Fixed-domain methods

derivation of (6.13) we find

I

Ti

{u.i.

(Ka<b)+H(U)

a<b}dX

dt

o 1

ax ax

at

1

8(0)

IT

a<b

= -

H(uo)<b(x,

0) dx -

K(gt)

ax (0, t) dt

+ 1 [{KU

a<b

-

K<b

au} dx -

H<b

dt] .

=.(1)-0

ax ax

A similar relation holds in R

z

.

By

adding these two relations, subtracting

(6.13)

and

assuming U is continuous across x = s(t) we obtain

1=8(t)

<b{[K

au/ax] dt+[H] dx}= 0,

(6.14)

where, in general,

[F]

denotes

the

jump

discontinuity in

the

function F

from x =

s(t)-O

to

x = s(t)+O. Since (6.14) holds for

an

arbitrary test

function

<b

and

an

arbitrary time, it follows

that

ds/dt =

-[K

aT/ax]/[H].

(6.15)

In

the

classical one-dimensional Stefan problem,

the

phase

change occurs

on

the

boundary x = s(t)

on

which U = Urn

and

from (6.1) we know

that

[H]

= pL. Thus (6.15)

is

seen

to

be

the

usual Stefan condition

[K

aT/axli =

-pL

ds/dt. (6.16)

We

have thus established, through

the

derivations

of

(6.13)

and

of

(6.16), that

both

the

classical solution

of

a

heat

conduction problem with

no

change

of

phase

and

the

classical solution

of

a one-dimensional Stefan

problem satisfy

the

criterion for a weak solution. Rubinstein (1971)

showed

that

the

classical two-phase Stefan

problem

(equations (1.7-9)

inclusive) is well

posed

with a unique solution provided

the

initial

condition

and

the

conditions

on

the

fixed boundaries

are

twice differenti-

able. Similar statements can

be

made

about

the

one-phase problem

(Sherman 1970; Tayler 1975) subject

to

further restrictions

on

certain

parameters. Furthermore, various authors have

proved

the

existence

and

uniqueness of

the

weak

solution (6.13) (Kamenomostskaja 1961; Oleinik,

1960) and

other

suitably defined weak solutions for multi-phase

and

higher-dimensional Stefan problems with derivative boundary conditions

(Friedman 1968). Rubinstein (1971) presented

the

uniqueness proof by

Kamenomostskaja

and

the

existence proof

by

Oleinik for a simple case.

A relevant survey

by

Niezg6dka (1983) and

more

references

are

discus-

sed

in §6.4.7.

From all

the

various studies

of

classical

and

weak

solutions we can

conclude

that

for any heat-conduction

problem

for which a classical

Enthalpy method

223

solution exists,

the

unique weak solution which exists will also

be

the

unique classical solution.

In

§6.2.2 it will be seen that certain finite-

difference schemes can

be

proved

to

converge

to

the unique weak

solution (Oleinik 1960) despite the discontinuities in

the

derivatives of u

across

the

phase-change boundary.

The

existence of

the

weak solution is

thus established through

the

proof of convergence.

The

discontinuities

make it difficult

to

prove directly

the

convergence of finite-difference

solutions

to

the classical solutions.

The

weak solution, therefore, has

important properties

both

theoretical

and

numerical. Indeed Rubinstein

(1979) claimed that it was only for

the

weak solution of

the

multi-phase

Stefan problem in several space variables that existence

had

been proved

and

saw this as a general position in

the

theory of partial differential

equations.

He

was critical of attempted proofs by Datzev (1970) and

Budak

and

Moskal (1970, 1971). NiezgOdka (1983) however, consi-

dered that

the

existence of classical solutions for some multi-dimensional

problems

had

been established by Meirmanov

(1981a,b-all

in Russian).

Considerable progress has certainly been made in theoretical studies of

one

and

two-phase multi-dimensional problems

by

formulating them as

variational inequalities (see §§6.4-6.4.7).

6.2.2.

An

explicit finite-difference scheme

As

a simple example, we describe an explicit finite-difference scheme

for

the

solution of

the

one-dimensional form of (6.9), where the variable

v is related

to

temperature,

u,

through (6.8). Standard modifications will

allow

the

scheme

to

be

adapted

to

include the form (6.11) with constant

or

variable

K,

for example,

and

also alternative conditions

on

the

fixed

boundaries.

The

derivative conditions

and

the

initial condition,

iJv/iJx=av,

iJv/iJx=O,

x=O,

x=t,

H(x,

0) = ho(x),

(6.17)

(6.18)

(6.19)

studied by Atthey (1974) provide a convenient introduction

to

his more

general problem with body heating discussed in §6.2.6. Also,

the

numeri-

cal scheme will require (6.10)

to

be

rewritten

v=H/d

b

H<O;

v=(H-Lp)/d

2

,

H>Lp;

v=O,

O~H~pL,

(6.20)

where, for convenience,

d;

=

pc;/

K;,

i = 1, 2.

The

simplest explicit replace-

ment of (6.9) is

(6.21)

224 Fixed-domain methods

where A = 8t/(8x)2,

8x

=

elM,

8t

=

T/N

for positive integers M,

N;

T>O

is

an arbitrary

upper

time boundary; and H;:' =

H(m

8x, n

8t)

etc., with

1-=:;m-=:;M-1,

O-=:;n-=:;N.

The

conditions (6.17-19) are taken as

v(j=vl'(2-a8x)/(2+a8x),

VM-l

=VM,

O~n-=:;M,

O-=:;m-=:;M.

(6.22)

(6.23)

The

first of (6.22) comes by approximating (6.17) as

(vl'-v(j)/8x

=

~a(vl'+v(j),

where

Vo

is

the

value

at

x=O.

Assuming

the

numerical solution has progressed

as

far

as

t = n 8t,

so

that values of v;:', H;:' are known for 0

-=:;

m

-=:;

M, we proceed

to

(n

+ 1)

8t

as follows:

(i)

compute H;:.+1 from (6.21) for

l-=:;m-=:;M-1.

(ii)

derive

v;:.+1

from

one

of the three relationships in (6.20) appro-

priate

to

the

value

of

H;:.+l.

(iii) compute

v(j+1

and

V;:'+l from (6.22).

Clearly this enthalpy scheme avoids

the

difficulties of front-tracking.

Instead, the position of

the

moving boundary has to

be

determined in

retrospect by inspection of

the

computed values of

Hand

v.

Since (6.21)

is a conservative integrated form of (6.9) for

the

element

8x

in

the

interval 8t,

the

acute problem of

the

infinite qerivative iJH/iJt is avoided.

Nevertheless,

the

numerical solution is likely

to

exhibit a stepwise be-

haviour and this may impair

the

accuracy with which

the

boundary can

be

located. Ways of improving this accuracy are considered in §6.2.5,

together with comments

on

implicit finite-difference schemes and ex-

tended applications of

the

enthalpy method.

Before leaving

the

simple explicit scheme (6.21) we discuss

the

con-

vergence of

the

numerical solution to

the

unique weak solution.

The

details of

the

outline of

the

convergence proof which follows are

to

be

found in papers

by

Kamenomostskaja (1961), Oleinik (1960), Atthey

(1972, 1974, 1975).

Piecewise-constant functions

f5

and H are defined by

f5(k,

t) =

v;:',

H(x,

t) = H;:',

(6.24)

for

(m

-1)

8x-=:;x

<m

8x, n 8t-=:;t-=:;(n+ 1) 8t, for each m, n in

the

numeri-

cal solution. Kamenomostskaja

(1961), having established

the

uniform

boundedness of

v;:', appJied relevant theorems of functional analysis

to

prove that there are sequences of

the

integers

M..

N

r

~

00

as r

~

00

such

that the corresponding subsequence of

the

sequence

{f5}

of functions v,

defined in (6.24), converges strongly

to

a function

v.

Also, the corres-

ponding subsequence of

{H}

converges weakly

to

a function H

as

r

~

00.

Enthalpy method

225

We

can now show

that

these functions v and H form a

weak

solution of

the

phase-change problem.

We

define, for any test function

<1>,

any choice

of integers M, N, and corresponding step-lengths

Sx,

St,

<1>;:'

=

({m

-n

Sx,

n St),

m=l,

...

,M,

n = 0,

...

,N,

(6.25)

<1>0=

<1>7(2-0:

Sx)/(2+0:

SX),}

n = 0,

...

,N,

~+1

=~,

together with associated piecewise-constant functions

(6.26)

cf,

=

cf,;:.,

cf,t

=

(;:'

-r:.;l)/St,

(6.27)

for (m

-1)

Sx';;;;x

<m

Sx,

n

St.;;;;t«n

+

1)

St.

We

can verify

that

cf"

<bx,.,

and

~

converge strongly

to

<1>,

ii<l>/ax

2

,

and

a<l>/at

respectively as

M,

N

~

00.

Now we proceed as

in

the

development of

the

weak solution (6.13)

and

multiply (6.21) by

<1>;:'

and

sum

over m and n

to

obtain

M

N-l

Hn+1

Hn

M

N-l

n 2 n + n

St

Sx

L L m - m

<1>;:'

=

St

Sx

L L

Vm+l

- v

m2

V

m

-l

<1>;:'.

(6.28)

m=l

n=O

St

m=l

n=O

(Sx)

On

summing

by

parts, inserting

the

initial and boundary conditions, and

finally replacing

the

sums

by

integrals we find

IT

r

ilcf,t

dx

dt

+ 1 T r

vcf,xx

dx

dt

= - r ho(x )

cf,

(x, 0) dx, (6.29)

where ho(x) is a piecewise-constant approximation

to

the

function ho(x)

in

(6.19).

In

the

limit as M

~

N

~

00

through

the

values

Mr

and

N..,

with

which subsequences

of

{v}

and

{if}

converging to v

and

H respectively

have

been

associated above, (6.29) becomes

IT

r{H

a<l>/at+v

a

2

<1>/ax'1)

dx

dt

= - r ho(x)

(x,

0) dx. (6.30)

Thus

the

functions v

and

H satisfy a weak solution.

The

final requirement is

to

show

that

the

solution

of

the

finite-

difference scheme (6.21) converges

to

the

weak solution as M

and

N

~oo

through any sequence of integers for which

the

relation A =

St/(Sx)2.;;;;

Ao<! is satisfied.

If

We

suppose this is not so

then

for some B

>0

and

for

some function,

I/J,

which is twice continuously differentiable, we

may find subsequences

of

the

sequences

{il}

and

{il} denoted by

{ti}

and

'"

{H}

for which

(6.31)

226 Fixed-domain methods

But

the discussion above about

the

properties

Qf

the

sequences

{v}

and

{il}

apply equally

to

the

sequences

{il}

and {H}.

Thus

we may find a

subsequence of

{il}

which converges strongly

to

a function

VI

and a

'"

subsequence of

{H}

which converges weakly

to

some function

HI'

Now

the

arguments used

to

establish (6.30) can

be

used

to

show that TI and

HI

form a weak solution.

But

we know the weak solution

to

be

unique

and

so we have

IIv

-

VIII

=

0,

IIH -

HIli

=

O.

Thus the above subsequences of

{il}

and {Ii} must converge

to

V

and

H, which contradicts (6.31).

We

have

therefore shown

that

the

solution of

the

finite-difference scheme (6.21)

converges

to

the unique weak solution (6.30).

Convergence proofs for

both

explicit

and

implicit finite-difference

schemes (see §6.2.4)

and

for discontinuous

and

smoothed enthalpy func-

tions (see §6.2.3) have been given by

other

authors including Ciavaldini

(1975), Jerome (1977), Schafer (1977), and Meyer (1973).

6.2.3. Alternative forms

of

the enthalpy function

One

definition of

the

enthalpy relation

H(u)

has

been

given in (6.1).

For

some substances, however, the phase change occurs within a wide

band

of temperatures

and

the

relationship between

H(u)

and tempera-

ture

u is known from experimental

data

or

standard physical tables. This

information can

then

be

used directly in the numerical solution (see, for

example, Albasiny 1956; Longworth 1975; Hodgkins

and

Waddington

1975).

Because of

the

difficulties and inaccuracies associated with the discon-

tinuous jump in

H(u)

at

the melting temperature

in

the

definition (6.1),

which is likely

to

occur for

pure

substances, various authors have prop-

posed

smoothing over a small temperature zone. With reference

to

the

definition (6.8) over

the

range - e

,.,.;;

v

,.,.;;

e where e > 0 is small Furzeland

(1977b) lists three examples of linear smoothing defined by

{

a2V-~'

v<-e,

H(v)=

a3v+L3,

-e"";;v"";;e,

alv+Lt.

v>e,

(6.32)

where the smoothing parameters

a;,

1..;,

i = 1, 2, 3,

are

chosen

by

Meyer

(1973) to

be

al

= d

l

, a2 = d

2

,

LI

= A

-(CI

+c2)e,

(6.33)

and by Elliott (1976)

to

be

al

=dt.

a2=d

2

,

a3=!{e(al+a2)+A}/e,

LI

=A,

L

2

=0,

L3=!{e(a2-al)+A},

(6.34)

Enthalpy method

227

where A = pL,

the

latent

heat

per

unit volume.

For

e = 0

the

smoothing

tends

to

a2 = d

2

,

a3

=

~(al

+ a:0,

L

2

=O,

L3=~A.

(6.35)

Friedman (1968) proved existence of a weak solution

to

the

Stefan

problem by showing that solutions for smoothed enthalpy converge

to

the

solution of

the

Stefan problem as e

~

o.

No

theoretical analysis has been made of

the

effect of

the

smoothing

parameters

on

the

accuracy of

the

solution.

The

influence of the width of

the

smoothing zone has been demonstrated numerically in selected cases.

Thus Bonacina

et

al.

(1973), whose findings were later confirmed by

Furzeland (1980), found their results for a one-dimensional, two-phase

Stefan problem

to

be

appreciably dependent

on

the

zone width, 2e, in

(6.32).

In

fact, Bonacina et

al.

(1973) replaced

the

enthalpy jump

at

the

phase-change by an equivalent heat capacity,

C(u),

i.e. density x specific

heat. Instead of an enthalpy

H(u)

defined by

f

u

{1,

u~Ur,

H(u)

=

C(u)

du

+

A'Il(U

-

Uf),

'Il(u) = 0

u,

,

u<Ur,

where

u,.

<

Ur

(freezing temperature) is arbitrarily chosen, they defined

C(u)

by

-

dH(u)

C(u)=~=C(u)+A8(u-Ur),

C(u)

=

{C

1

(u),

u<Ur,

C2(u), u >

Ur,

where

8(u

-

Uf)

is the Dirac function.

The

two-phase Stefan problem is

then

expressed by

the

single equation

-

au

a {

au}

C(u)-=-

K(u)-

,

at

ax ax

Bonacina et

al.

(1973) then smoothed

the

Dirac function in

the

definition

of

C(u)

and

also in

K(u)

and adapted a linear form for both

C(u)

and

K(u)

over

the

range

Ur-e

~u~uf+e.

In

a similar comparison, Shamsundar (1978) found

the

heat-flux-time

curve

to

be

relatively smooth

but

displaced from the

true

curVe for a large

value of

e;

as e was reduced

the

smoothed-enthalpy calculated curve

approached

the

correct result

but

developed a stepwise wavy character.

Elliott (1976) quoted an example in which an optimum choice of

e can

improve

the

numerical solution as much as a reduction in mesh size.

Furzeland

(1977b) quotes results

but

points

out

that

the

optimum e has

228 Fixed-domain methods

to

be

determined by experiment. Meyer (1973)

on

the

other

hand

claims

his calculated values for

the

temperature record at one point in

the

region

for a two-dimensional melting problem around a square duct are inde-

pendent of

B in

the

range

1O-6~B

~0.5.

Incidentally, Meyer (1973)

mentions a precipitation problem (Cannon and Hill 1970) in which

the

enthalpy actually is continuous and piecewise linear.

Budak

et

al.

(1965) and Moiseynko and Samarskii (1965) used higher-

order smoothings in which, for example,

the

H(u)

function

near

the

melting temperature is approximated by a parabola which ensures con-

tinuity

of

the

first derivative

of

the

approximating smoothed enthalpy

function

as

well as of

the

function itself. Brauner et al. (1983) discuss

the

smoothing function

H(u)=!Lp{l+u/('7I+lul)},

where '71>0 is a small

parameter.

6.2.4. Other numerical schemes and multi-dimensional problems

The

discontinuous form of

the

enthalpy function expressed

as

u =

u(H)

was used by

Atthey

(1974) together with

the

one-dimensional explicit

form (6.21) though

he

did

not

use

the

transformation (6.8).

The

scheme

was extended

to

two-dimensional problems by Crowley (1978).

She solved two problems

of

the

inward solidification of a square cylinder

of liquid initially at its freezing temperature for two different surface

conditions. This problem was described in §5.4.2(i).

In

the

first case,

the

surface temperature is lowered

at

a constant rate, corresponding

to

the

conditions under which Saitoh (1976) carried

out

his experiments. Crow-

ley's (1978) enthalpy calculations using an explicit finite-difference

scheme agree well with Saitoh's experimental results and with a perturba-

tion solution in inverse powers

of

the

latent heat,

L,

until

the

freezing

front has moved about halfway

to

the

centre.

In

the second case,

the

surface temperature is dropped discontinuously

at

the

initial instant.

The

solution obtained by

Rathjen

and

Jiji (1971) for

solidification in

an

infinite corner provides a good approximation while

the

front remains parallel

to

the

sides of

the

square cylinder away from its

corners. Crowley also compared

her

enthalpy results graphically with

those of Allan and Severn (1962), Lazaridis (1970),

Crank

and

Gupta

(1975). .

Implicit finite-difference schemes with discontinuous enthalpy functions

were described

by

Furzeland (1974, 1977b), Federenko (1975), Long-

worth (1975), and Shamsundar and Sparrow (1975). Wood, Ritchie, and

Bell (1981) used a hopscotch finite-difference scheme (Gourlay, 1970).

Ciavaldini (1975) 'used explicit and implicit finite-element schemes

to

solve a discretized weak form

of

an enthalpy formulation. Hodgkins and

Waddington (1975) also introduced finite elements.

Li

(1983) adopted

the

finite-element discretization used by Bonnerot

Enthalpy method

229

and

Jamet

(1979), shown in Fig. 4.5,

but

applied it

to

a

weak

enthalpy

formulation. This is a front-tracking method which avoids explicit treat-

ment

of

the

jump

condition

and

so overcomes

the

difficulties

of

an

accurate, a posteriori location

of

the

moving boundary discussed in

§6.2.5.

The

method

is applied

to

a one-dimensional, two-phase problem,

a two-dimensional, one-phase problem, the spot-welding problem

(§6.2.6),

and

to

alloy solidification (§6.2.7(ii)).

Longworth (1975) considered

an

equation of

the

form

...

aHtat+u'

VH

-

VZcf>

= Q(x, t), (6.36)

which describes a Stefan melting problem with a

heat

source Q(x, t)

moving with a velocity u,

and

cf>

is

the

Kirchhoff variable defined

in

eqn

(6.8).

The

equation refers

to

a frame

of

reference fixed in

the

heat

source,

and

Longworth integrated over individual, elementary, volume cells to

obtain for

the

ith cell

of

volume

V;

and surface area

SI

aHJat = Q

i

-J.

Hu'

dS

+J.

Vcf>·

dS,

s. s.

where Q

i

=

(ltV;)

Iv, Q dV.

The

surface integrals are expressed

in

terms

of

H;,

cf>i

and

~,

cf>j

for cells j adjacent

to

cell

i,

but

Longworth expresses

all

the

cf>s

in terms

of

the

corresponding

Hs

through

the

enthalpy

relationship

and

thus arrives

at

a

set

of

non-linear finite-difference equa-

tions

of

the

form

aB/at = F(H),

(6.37)

where H is a vector with components H; and F is a non-linear function

whose

ith

component

Fi

is a function of H; and

the

values

~

for

the

cells

j adjacent

to

cell

i.

An implicit Crank-Nicolson scheme is evaluated

iteratively.

Smoothed enthalpy functions were first used in multi-dimensional

problems

by

Budak

et

aZ.

(1965) and Moiseynko and Samarskii (1965),

with locally one-dimensional finite-difference methods. Two-dimensional,

implicit finite-difference and finite-element schemes were suggested by

Couch et

aZ.

(1970), Elliott (1976), and Meyer (1973, 1975, 1976,

1978c). Meyer had a particular interest in

the

efficient numerical solution

of

the

non-linear equations which

the

implicit formulations require.

Bonacina et

aZ.

(1973), Fisher

and

Medland (1974), and Comini et

aZ.

(1974) used three-level schemes.

The

latter authors aimed particularly

to

include simultaneous temperature dependence

of

thermal conductivity,

heat

capacity,

rate

of

internal

heat

generation, and surface

heat

transfer

coefficients.

Latent

heat

effects were approximated by a large heat

capacity over a small temperature interval enclosing

the

melting tempera-

ture. An improved procedure was described by Morgan et

aZ.

(1978).

230 Fixed-domain methods

For

the

implicit schemes,

there

arises

the

complication of having

to

decide in which

of

the

three enthalpy intervals a given point lies

at

the

next time level in

order

to

use

the

correct relation from (6.1)

or

(6.32),

for example.

The

authors cited above determined which interval was

appropriate by a trial and

error

process, and this was time consuming

especially as it was

an

additional

part

of

an

iterative scheme

of

solution.

Furzeland (1977b) described a way

of

avoiding

the

trial

and

error

process

that

he

attributed

to

Shamsundar

and

Sparro,¥ (1975).

A fully implicit, finite-difference scheme approximating

the

two-

dimensional case

of

(6.9)

is

H

n+1+4

n+1_Hn

+ {

n+1

+

n+1

+

n+1

+

n+1}

;'i rvi,i -

iJ

r

Vi.J-1

Vi-1.i

Vi+1.i

V;'i+1'

(6.38)

where

h=8x=8y,

r=8tlh

2

,

vi;J=v(ih,jh,

n8t),

HfJ=H(vi;i)'

The

Gauss-

Seidel iterative scheme, working along

the

rows j

and

up

the

columns

i,

can

be

written

where

the

notation n + 1 has been dropped for

the

unknowns and

b

- {

(s+1)

+

(s+1)

+

(s)

+

(s)

}

;'i - r Vi.i-1 Vi-l.i

Vi+1.i

vi.i+1,

so

that

b

ii

is known

at

the

iteration (s + 1).

(6.39)

(6.40)

With reference

to

the

temperature-enthalpy relations (6.20), H

<0

implies v < 0 and so

bi.J

<

O.

Conversely,

if

b

i

•

i

< 0,

then

v is given by

the

first of (6.20) and using this

in

(6.39) gives

H~s.+1)

= b

i

•

i

'J

1

+4rld

1

if

Similarly

H~s.+1)

= b

i

.i

+4rAld

2

'.J

1+4rld

2

'

bi.i>A,

(6.41)

(6.42)

(6.43)

where A

=!.p.

Thus,

the

known value

of

b

i

•

i

at

any

point

in

the

computa-

tional scheme determines which expression in (6.20) is

to

be

used.

The

argument is only valid as long as

the

dependence

on

K has

been

removed,

e.g. by

the

transformation (6.8),

and

provided

both

Hand

U are negative

in

one

phase

and

positive in

the

other. This second requirement can easily

be

satisfied for cases where

the

melting temperature

is

Urn

1=

0 by adopting

new variables

u'

= U - Urn,

H'

= H -

PC2Urn'

Finite elements can

be

used

instead

of

finite differences

and

Furzeland (1977b) includes a body-

heating term Q(x, y) which, by

the

simple modification

of

including Q

on

the

left-hand side

of

(6.39),

can

be

extended

to

Q(u).

Enthalpy method

231

The

same arguments

can

be

applied (Furzeland, 1977

b)

for the

smoothed enthalpy function (6.32). H

v<-e,

then

H<-(a2e+L~

and

so H +

4",

<

-{

e(

a2 + 4r) + L:J. This corresponds

to

the

test (6.41) for

L2~0

and, in fact, (6.41), (6.42),

and

(6.43)

are

to

be

replaced by

(.+1)

_ bi•

i

+

L2

b

Vi,i

- a2+

4r

'

bi,i<-

2,

(6.44)

(.+1)_

bi,i-

L

3 b b b

v··

- -

2:S:;

iJ:S:;

1>

',J a3+

4r

'

(6.45)

(0+1)

_

b;,i

-

Ll

b b

Vi,i

-

al+

4r

'

i,i>

1>

(6.46)

where b

1

=e(al+4r)+L1> b

2

=e(a2+4r)+Lz.

Furzeland (1977 b) applies this iteration -avoiding modification

to

obtain

fully implicit, finite-difference solutions

of

the

constructed problem in two

space dimensions (x, y):

c;

au/at

=Kt

V

2

u+ Q"

i = 1, 2,

with

Q

i

= Qi(t)

=2(C;e-

2t

-2Kt),

and Q

1

is used

if

x

2

+y2>e-

2t

, Q

2

if

x

2

+y2<e-

2t

.

The

initial conditions

are

t=O,

where s(x,

y,

t)

= 0 is

the

moving boundary

on

which u = 0 and

the

usual

Stefan jump condition holds.

The

region

x~O, y~O,

x2+y2:S:;1 is consi-

dered

with

the

boundary conditions

au/ax

=0,

x=O;

au/ay

=0,

y=O.

The

solution

of

this problem is u = x

2

+y2_e-

2t

with

the

moving bound-

ary given by s(x, y,

t)=x

2

+y2_e

2t

=0.

Numerical solutions obtained for

both

discontinuous and smoothed enthalpy with

the

parameters (6.33)

and (6.34) are reproduced in Table 6.1.

Furzeland (1980) outlined

the

basis of a general-purpose program by

Elliott

and

Furzeland (Furzeland, 1979b) for

the

finite-difference solution

of

the

non-linear, one-dimensional problem

i=1,2,

(6.47)

(6.48)

Crank J. Free and Moving Boundary Problems

Подождите немного. Документ загружается.