Crank J. Free and Moving Boundary Problems

272

Fixed-domain methods

equations (6.160) to (6.164)

and

leads to the equations

ue~

-V,

(KVu)-t)=O

in

D,

U~O,

u(x,O)=O in D,

au '

--V,

(KVu)-t~O

in D

at

u=g(t)

on r

o

,

where

g(t) = f g(T) dT

u=O

on r

c

,

for

u=Kau/an

onr

F

,

and

&(0)=0.

(6.167)

So far,

the

initial state has been

tiken

to

be

water

at

the

freezing

temperature zero everywhere i.e.

6(x, 0) = 0 in D. Ichikawa and Kikuchi

(1979) use a transformation due

to

Friedman and Kinderlehrer (1975) to

extend their treatment to

the

initial state which

is

partially ice imd

partially water

at

zero temperature. Taking 0

0

to denote the ice domain

at

the

initial stage, the freezing index u

is

redefined as

u(x, t) =

it

6(x,

T)

dT

in n-Oo,

sex)

= f 6 (x,

T)

dT in

00,

.

=0

in

D-O.

The

partial differential equation for u then becomes

au/at=v·

(KVu)+h,

(6.168)

(6.169)

where

h = t

in

0-0

0

and h = 6

0

(x)

in

00,

where 6

0

is

the initial

temperature in

00.

Ichikawa and Kikuchi (1979) further point

out

that

if

r F

is

part of

the

boundary of the domain 0

0

just defined,

the

more general boundary

condition '

K a6/an =

0:6+13

on

r

F

can

be

considered, and that on integration with respect to time from t = 0

the corresponding conditioll for

u

is

•

K aujan =

au

+

I3t

on r

F'

If

r F

is

loc?ted in the boundary of

00

and

0,

the

generalized relation

cannot

be

used.

By

multiJ-lying

the

expression

au/at-V'

(KVu)-t

Variational inequalities

273

by (v - u) and integrating by parts we obtain a generalized inequality,

.which incorporates the conditions

on

r

F,

(Ur,

v-u)o

+(KVu,

V(v-

u»o-(t,

v-

u)o =

(Ur-V·

(KVv)-t,

v)o

-(au,

v -

-u)rF;;'-(au,

v -

u)r

F

.

(6.170)

Here

the test function v ~ 0 is used

in

D for all time t in 0 ~ t ~ T with

v = 0

on

r c, v = g

on

r

G,

where T is a positive real munber indicating the

total time of

the

problem. Also, the inequality

(Ur-V·

(KVu)-t,

v)o;;'O

has been used and

v~O

so that

the

modified form of (6.170),

(Ur,

v -

U)D

+(KVu,

V(v -

U)D

+ (au, v - u)rF;;.(t, v - u)o, (6.171)

is obtained. Ichikawa

and

Kikuchi (1979) quote theoretical work by Lions

(1976) and Ichikawa (1978) which establishes the existence of a unique

solution

of

(6.171) for functions v which, together with their first deriva-

tives, are square summable

on

D and satisfy

the

boundary conditions as

above.

Ichikawa and Kikuchi (1979) discretized in space using finite-elements

and

a generalized implicit scheme in time. They wrote (6.171) as

(Ur,

v -

u)+B(u,

v - u);;.(f, v - u)

for v as defined above an'd where

B(u,

v) =

(KVu,

VV)o + (au,

v)r

F

,

(u,v)=(u,v)o,

(f,v)=(t,v)o'

(6.172)

The

derivation corresponding

to

that

of (6.147) shows

the

discretized

form

of

(6.172) to

be

N KijUj+1(Vi -

u?+1)

;;'/;(tlj -

U?+1),

(6.173)

where

u = L

U;<Pi

=

u;'<Pb

f =

t<Pb

V =

WA,

i=l

and

The

projectional

SOR

algorithm becomes

u;:J1(k+!) = (1-w)u;:.T

1

(k)+w(

_if

K;.ju;:.?(k+1)

j=1

~

A

+1

A)

A

-

t...

~.jU;:.j

(k)

+t

/~i

j=i+l

u;:t1(k+ 1) =min(O, u;:t

1

(k+!»,

(6.174)

274 Fixed-domain methods

y

10

G

FED

" C

I

\

'-

B

O~----

__________

~A~_x

10

FIG. 6.10. Two-dimensional problem with coupling freezing fronts

where u/:,tl(k) denotes the kth iteration of u/:,t

l

. Ichikawa and Kikuchi

(1979) solved

the

variational inequality (6.173) for

the

one-dimensional

single-phase Stefan problem using linear finite elements with different'

mesh sizes and found the orders of

error

to

be exactly

the

step size h in

the

HI

norm and

h1-6

in the

L2

norm in agreement with their theoretical

estimates. They also quoted results obtained by several schemes for

solving

the

inequality for a model, one-dimensional, steady-state problem

with a known exact solution and concluded that

the

Langrange multiplier

method was unsatisfactory and that

the

projected

SOR

method

~as

more

effective than either

the

fixed point

SOR

or

the penalty method.

Numerical experiments on a two-dimensional problem illustrated iii

Fig. 6.10, and selected because the two freezing fronts become coupled

TABLE

6.12

Coordinates

of

the freezing front on the axes

LCO'

JMMh

8x

==

1/20 8x = 1/80 8u = 0.2 8u = 0.1

0.05

0.776

0.770

0.775 0.775

0.10 0.680

0.673

0.675

0.676

0.20 0.539 0.537

0.535 0.536

0.30 0.431 0.426

0.422 0.420

0.40 0.323 0.321 0.314

0.308

0.50 0.208 0.209

0.201

0.188

a

LC

=

linear

complementarity method (Elliott and

Ockendon 1982);

8t=i(8x)2.

b IMM = Isotherm migration method (Crank and

Crowley 1978);

8t

5

0.OOOl.

Variational inequalities

275

TABLE

6.13

Coordinates

of

freezing front on diagonal

y=x

Le"

IMM"

8x=1/20 8x=1/80

8u=0.2

8u=0.1

0.05 0.728

0.735 0.732 0.732

0.10 0.625

0.625 0.623 0.619

0.20 0.470

0.470 0.471 0.463

0.30 0.350 0.351

0.352 0.342

0.40 0.247

0.249 0.250 0.237

0.50 0.152

0.154 0.153 0.139

• See footnote

to

Table 6.12.

after a finite time, reveal an optimum relaxation factor of w = 1 in the

projected

SOR

method used. Lumping of the mass matrix gives smooth

freezing fronts in

the

first few time steps, whereas singular fronts are

observed in the absence of lumping.

The

two solutions agree

at

later

times. '

. Elliott (1976) also quoted numerical results, obtained by solving the

variational-inequality formulation of one-phase Stefan problems, which

agree satisfactorily with

the

analytical solution. His solutions for the

melting of a long cylinder of square cross-section, quoted by Elliott and

Ockendon (1982), are compared with those obtained by Crank and

Crowley (1978) using

the

isotherm migration method described in

§5.4.2(iii) in Tables 6.12, 13.

6.4.6. Two-phase Stefan problem

As an example of the extension of

the

freezing-index transformation to

a two-phase Stefan problem and

of

differential and variational formula-

tions in multi-space dimensions, we refer

to

a paper by Kikuchi and

Ichikawa (1979). Their main concern was to present numerical schemes

of solution based

on

theoretical considerations due

to

Duvaut (1975),

Fremond (1974), and Lions (1976).

Consider a multi-dimensional domain

D to be divided into a solid part,

Db

and a liquid part, D

2

,

with temperature u(x,

t)<O

in D1 and u(x,

t»

o in D

2

•

Here

x

is

a space vector of position with one, two,

or

three

components.

On

the phase-change interface f o(t) we have u(x, t) =

O.

We

take

a problem in which there may

be

different conditions on three

separate parts of

the

boundary f of D. The temperature u(x, t)

is

prescribed

on

the boundary f

1

; there

is

no heat

flux

on f

2

; on

f3

the heat

flux

is

taken

to

be

proportional

to

the

temperature.

The

problem can

be

defined

276

Fixed-domain methods

by

the

equations and conditions

C

1

iJU/iJt

= V .

(K

1

Vu) in

Db}

C

2

iJU/iJt

= V .

(K

2

Vu) in D

2

,

u(x,

t) = g(x, t)

on

f1'

}

iJju(x,

t)=O

on

f2

bounding D

..

iJju(x, t) =

-pju(x,

t)+qj(x,

t)

on

f3

bounding D

..

u(x,

t)=O

on

fo, }

[KVu]i'

Vl+L

=0

on

fo,

u(x,

0)

= uo(x) in D.

(6.175)

(6.176)

(6.177)

(6.178)

Here

G and K

..

i = 1, 2, are

the

heat capacity

and

conductivity of

the

ith

phase,

g,

P

..

qj are" given functions,

Uo

is

the

initial temperature in

the

domain

D,

iJiu

is

the

temperature gradient normal

to

f,

and

[KVu]i=

KiVu)2-

K

1

(Vu)1'

The

function

t=

lex)

indicates

the

position

of

r

o

, i.e.

x = set) as in §6.4.5.

For

a two-phase problem Duvaut (1975) introduced

the

transformation

w(x,

t) = r

K;u(x,

T)

dT.

(6.179)

For

u(x,

t) continuous

we

haveK;u(x,

t)=iJw/iJt=

w(x,

t)

and

so w(x,

t)<O

in D

1

(t)

and

w(x,

t»O

in D

2

(t).

On

integrating (6.175) with respect

to

time

and

remembering

that

a point x may

be

in D

1

(0)

initially

but

in

D

2

(t)

at

time t we obtain relationships such as

<;

V . Vw(x, t) =

-C

1

u(x,

0)+

C

1

u(x, t), x in D

1

(0) and D

1

(t),

=-C

1

u(x,

0)+

C

2

u(x,

t)

+ L, x in D

1

(0).and

Dit),

and

two corresponding relations where x is initially in D

2

(0)

but

in D1 (t)

or

D

2

(t)

at

time

t.

Using (6.178)

the

equation which holds over

the

whole

domain D can

be

written

where

iJw

(GlK;)

- (x,

t)-V'

Vw(x,

t)=

Cjuo(x)+ EjjL

iJt

i = 1,

x'

in D

1

(t),

j = 1, x in D

1

(0),

.

E12

= 1,

E21

=-1,

i

,:,,2,

x in

Dit),

j = 2, x in D

2

(0),

Ell

=

~22=0

.

(6.180)

(6.181)

(6.182)

(6.183)

Now

iJw

=

J~K;

Vu . n

dT

=

J~

iJju

dT,

where n is

the

outward normal unit

Variational inequalities

277

vector

on

f,

so

that

boundary conditions (6.176) become

w(x, t) = g(x, t)

on

f

1,

aw(x,

t)=O

onf:z,

aw(x, t) =

-(VdKi)w(x,

t)+q(x,

t) .

on

f3'

(6.184)

(6.185)

(6.186)

where

i = 1 for x in

Dl(t)

and i = 2 for x in D

2

(t)

and

we have

put

g(x, t) = r Kig(x,,.) d,.

with

·e

if

g(X,T)<O'}

z=

if

g(x,

,.»0,

2

q(x,

t)=

f q;(x,,.) d,.

with

·e

if

x in Dl('T)

z=

2

if

x in

Di'T)

(6.187)

The

two-phase Stefan problem is now expressed over

the

whole Jixed

domain

D by

the

equation (6.180) with the conditions (6.184-6).

The

interface conditions (6.177) and

the

initial condition w(x,

0)

= 0 have

been included.

Under

the

assumption

that

w(x, t) is twice differentiable in D for all t

in

O~t~T,

Kikuchi and Ichikawa write (6.180) in

the

variational form

(d;wj,

v)+a;(w,

v) =

(Bi~'

v)+

Fj(v),

for all t in

O~t~T,

where

Fj(v) = L

quov

dx +

1.

qv ds,

a;(u, v) = L

VU'

Vv

dx+

1.

ei

uv

ds,

(6.188)

the

test functions v

and

their first derivatives are square summable in D

and

v(t)=O

on

fl'

Kikuchi and Ichikawa (1979) obtained finite-difference solutions based

on

a discretized form of (6.180) which in two space dimensions is

di,(:Iw",,(:I

+

(4l6?)wa,(:I

= (1/a

2

)(w",_1,(:I

+ W",+1,(:I +

Wa,(:I-l

+ Wa,(:I+l)

+

Cj'(:I(uo)a,(:I

+ Bij,(:IL (6.189)

where

d;

= CJKb

(0:,

(3)

denotes

the

point

(0:

ax,

(3

ay)

on

the

rectangular

net

covering

the

domain D and

Il.

=

ax

=

ay

is adopted for simplicity,

wa,(:I denotes

the

value of W at

the

point

(0:,

(3),

etc.

The

term

di'(:Iwa,(:I

is

approximated

by

di,(:Iw"',1l (n

at)

= di,(:I(n at){W:'(:I-

W:;l}/at,

where

w:'(:1

=

w;"(:I(n

at)

and

at

is

the

given time interval, and either

278 Fixed-domain methods

df,13(n

dt)

or

df,13(n

-1)

dt

is

used.

The

solution proceeds in time steps

dt

by a modified

SOR

scheme

as

follows, where m denotes

the

mth

stage of

the

iteration and

df,13(n

dt)

is

used (a simple modification

is

needed

if

di"13(n

-1)

dt

is

used instead); w

is

the

iteration parameter and 0 < w < 2.

(i)

Put

w:,13(O)

=

W:~l.

(ii)

R..,13(m)

= HW:-l,13(m) + W:+l,13(m-1)+

W:,13-1(m)+

w:,13+1(m

-1)

+ d

2

Cj,13(uo)}

N

" (

)-

",13(

)L

di"13(n,

m){

"(

-1)-

n-l}

..

,13

m -

Bi,i

n,

m

dt

W

...

13

m W

..

,13,

w:,13(m)

=

(1-W)W:,13(m-1)+w{R:.

13

(m)+

~2

N:,i

m

)}.

(iii) Repeat

till

I

IW:,13(m)

- w:'im

-1)IJI

IW:,13(m)1

<

B,

where B is

the

tolerance of convergence of

the

iterations.

The

indices

i,

j are given by

i=1,

w:,13(m)<O,

j=1,

(UO)

...

13<O,

i=2,

w",13(m)

>0,

j=2,

(uo)

...

/3>O,

where

w",13(m)=(1/dt){w:,13(m-1)-w:~1}.

Kikuchi and Ichikawa (1979)

argue that the erm

Bij,13(n,

m)L

varies like a step function, whereas the

other terms in

the

iterative scheme change moderately and therefore in

order

to

secure convergence the relationship

IR:,13(m)-i

d

2

di"/3(n,

m)w",13(m)1

»Ii

d

2

Bij13(n,

m)LI

(6.19<,»

must be satisfied.

Two problems are solved in one space dimension.

In

the

first,

the

initial

temperature

is

given by

uo(x)=(x-~

for

x~!,

uo(x)=O

for

~<x<1,

and the boundary conditions are

g(O,

t) =

-!,

g(1, t) =

t,

so that for t < 1,

&(0, t) =

-K

1

t/2,

&(1,

t) = K2t2/2. Numerical experiments support

the

con-

vergence condition

2Kl

dt;:;;:

Ld,

which

is

an approximation based

on

(6.

190a).

It

is

deduced by noting that for the one-dimensional case

(6.190) becomes, for

n = 1,

10.5(w~_1(m)

+

w~+l(m

-1»-0.25

d

2

df(1,

m)w!(m)\»

0.25L d

2

•

(6. 190a)

and then concentrating attention on the point

x = a d

=!

and assuming M

small enough for

w!

=

w~

=

Uo

..

approximately. Since

Uo..

= 0

at

x

=!,

UO,

..

-l

=

-d

at x

=!-d,

uo,

..

+1

= 0

at

x

=!+d,

and recalling that

L

~t

wl-l

= K1u

..

-l

d'T

=

Kl

dtUO,

..

-l

=

-K

1

dt

d

we arrive

at

0.5Kl

d'

At»0.25

d

2

L

as

the convergence condition.

Variational inequalities

279

The

second problem, also

one

dimensional, includes

heat

parameters

appropriate

to

a silty soil with

20

per

cent

moisture content; initially

the

soil is unfrozen with temperature uo(x) = ax

2

, where a = g(b,

0)/b

2

,

x is

the

depth

through

the

soil,

and

g(b,O)

the

initial temperature

at

the

full

depth

b.

Boundary conditions

for

g(O,

t), g(b, t) are prescribed numeri-

cally. Numerical convergence studies support the approximate criterion

\!K1 a

2

g(0,

.at)\»

a

2

L, derived as in

the

earlier Problem.

Graphic

solutions were presented by Kikuchi and Ichikawa (1979) for a

two-dimensional problem defined in

the

square region

O~x~O.4,

O~y~

0.4, with C

1

= 0.5, C

2

=

K1

= K2 = L = 1.0, and zero initial temperature.

Boundary conditions

are

g(0,y)=(0.8-y)(2t-1),

g(x,O.4)=

(0.8 - x)(2t

-1),

g(O.4, y) =

O.s.Jy(l-

3t), g(x, 0) =

0.5.Jx(1-

3t).

In

this

example,

the

term

8ij,13(n,

m)L

is replaced by

the

expression

{h1

(e",-l,13

+ 8",+1,13 + 8",,13-

1

+ 8",,13+

1

) + h

28

a.13}L/(4h

1

+ hz),

where

8 =

8ij

for convenience of writing, and

the

influence

on

the

results

is examined for different

hi>

h

2

•

Also

some results

are

obtained with a

linearly smoothed phase boundary,

and

there is a mushy region.

Finally, Kikuchi and Ichikawa

(1979) used linear finite elements

to.

solve

the

variational equation (6.188) for a two-dimensional problem

defined with reference

to

Fig. 6.11 by

the

conditions:

(

)

_

()_{-1.0,

O<t~O

4 and

0.5<t~0.6}

u x, t - g t -

05

04

",;:::05

on

fl'

. , . <t...., .

On

other

parts

of

the

boundary

the

normal derivative of u(x, t) is' zero.

The

heat

parameters

are

K1

=K2=

C

1

=L=

1.0, C

2

=0.5.

The

Cryer

algorithm is used with

w =

1.4

as

the

over-relaXation parameter.

The

shape

of

the

frozen area just

'before

it vanishes depends markedly

on

the

width

of

the

smoothed, phase-change region

but

does

not

seem

to

effect

the

temperature

distribution in subsequent time steps.

Kikuchi

and

Ichikawa (1979) did

not

explicitly include

an

enthalpy

r

2

0.4

x

Fw.

6.11. Twp-dimensional two-phase problem

280

Fixed-domain methods

function in

either

of their papers.

Duvaut

(1975)

and

Elliott (1976)

derived a variational inequality for

the

two-phase Stefan problem which

incorporates

the

enthalpy function and its integral with respect

to

temper-

ature.

If

heat

capacity C

and

conductivity K

are

taken

to

be

piecewise

constants, a convenient definition of enthalpy

H(z)

is

H(Z)=O'.lZ,

z<O;

H(z)

=

0'.2Z

+L,

z>O;

O~H(z)~L,

z = 0, (6.191)

where

0'.;

=

G/Kb

i = 1, 2. Consider a domain D divided into two

parts

Dh

D2 by

the

melting surface

at

t

==

lex)

on

which

the

usual Stefan

condition holds.

The

problem

is

O'.i

au/at

==

VZu in Di(t), i = 1,

,Z,}

.

u = g

on

the

boundary aD,

u=uo<O,

H(uo) = H

o

,

t=O.

(6.192)

In

this formulation

the

variable u is related

to

temperature T by

u

=.C

K(6)

d6,

(6.193)

which allows for likely discontinuities in

K(T),

and

the

freezing index

becomes

w(x, t) = f u(x,

'T)

d'T

(6.194)

;;

in place

of

(6.179). Proceedipg as in

the

derivation of (6.180) above

we

see

that

t<l(x),

t>l(x),

which can

be

combined into

the

single equation

V

2

w +

0'.1

Uo

= H(aw/at) =

H(

w,)

in place of (6.180).

(6.195)

Because

H(z)

is discontinuous

at

z = 0

it

is useful

to

introduce its

integral

4>(z)

which

is

the

continuous function

4>(z)

=

rH(Z')dZ'=!0'.1(Z-?+~0'.2(z+)2+LZ+

(6.196)

where

z-=min(O,z)

and

z+=max(O,z).

It

follows readily

that

4>(Zl)-4>(Z2);;;,:H(Z2)(Z1-

z~,

(6.197)

Variational inequalities

281

for all real

Zl

and provided H(Z2) satisfies (6.191). Combining (6.195)

. and (6.197) yields

(6.198)

The

desired variational inequality follows by substituting Z =

v(x)

and

integrating (6.198).

It

is

(6.199)

where

cI>(v)

= J

4>

(v

(x»

dx

for all admissible

v.

Duvaut (1975) proved that there exists a unique

solution of (6.199) for which

w

:0,

t=O.

6.4.7. Mathematical

results

Some references have been made in this chapter

to

investigations of the

mathematical properties of classical, weak, and variational solutions of

generalized Stefan-like problems.

It

is

nQt

intended in the present section

to

discuss these aspects in detail but simply

to

refer interested readers to

some recenj surveys and

the

numerous references contained in them. One

up-to-date survey by Niezg6dka (1983) concentrates on recent and

lesser known results and extends earlier surveys by Friedman (1979) and

Primicerio

(1981a) with reference

to

the bibliography by Wilson, Sol-

omon, and

Trent

(1979). A brief historical introduction is given by

Pawlow (1981).

Niezg6dka (1983) listed 16 applications of practical origin with source

references which introduce various generalized features and have moti-

vated research into

the

evaluation of solutions and their theoretical

properties. His survey of mathematical results for generalized statements

of Stefan-like problems covers one-dimensional and multi-dimensional

problems, both single and multi-phase.

The

concept of the weak solution introduced by Kamenomostskaja

(1961) and extended by Oleinik (1960) was based

on

the Kirchhoff

transformation

and

the

enthalpy function. Kamenomostskaja considered a

two-phase multi-dimensional problem and proved the existence and

uniqueness of the weak solution by considering

the

convergence of

explicit finite-difference approximations

as

in §6.2.2 above. Oleinik

(1960) used a smoothing technique

to

prove existence and uniqueness for

a more general one-dimensional Stefan problem.

The

work of both these

authors was refined and extended by Friedman (1968).

The

more general

problem with parameters dependent on space, time, and temperature,

and with a derivative condition on the fixed boundary, was investigated

Crank J. Free and Moving Boundary Problems

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