Crank J. Free and Moving Boundary Problems

182 Front-tracking methods

At

the

nth

time level we can approximate

the

system (4.41-46) by

(

'F

')'_

U;

-Vj(X)

D

l\iUi

C;

l)t

.c';,

ul(b

1

)

= {3uHb

1

)+

al>

u~(lh)

=

{32U2(b~

+ a2,

i=

1, 2.

(4.47)

k1uHs)-

k2U~(S)+

A(S

-l)~-l)

=

IL3

where the prime indicates

d/dx,

l)t

=

t,.

-

t,.-I>

and where all parameters

are evaluated

at

t =

t,..

The

functions

Vi

denote

U;

at

the previous time

level

or

their linear extensions beyond their domain of definition in

the

neighbourhood of

the

moving interface. Thus

if

the

boundary

is

moving in

the

direction of x increasing, so

that

s(t,.»S(t,.-l), we use

Vl(X) =

Ul(s,,-l)+(X

-s,,-l)uHs,,-l),

x ;;;'s,,-l' (4.48)

A

full description of

the

method of lines applied

to

the discretized

system (4.47) was given by Meyer

(1978c). We require to find

the

positions of

the

moving boundary

at

successive times t =

t,.

and

to

solve

the sequence of second-order, ordinary differential equations.

The

dis-

cretization of

the

time derivative does not depend on

the

linear structure

of the equations (4.41-6), and

if

these equations were non-linear

the

usual 'shooting technique'

could'

be

used

to

solve them. For linear

equations and conditions

on

the

fixed boundaries, Meyer (1973, 1978c)

developed a convenient method of solution by introducing

the

well-

known Riccati transformations,

to

be

used

at

t =

t,.

;

Ul(X) =

U(x)<h(x)+w(x),

}

q,2(X)

= R(X)U2(X)+

z(x),

«>;

= k;(x,

t,.)u:,

i = 1, 2.

We obtain

U,

w,

R, and z by solving

the

initial value problems

dR

C2

R2

-=---

dx

l)t

k2'

dz

= R

z-

C2V2+F

dx

k2

l)t

2,

where all parameters take their values

at

t =

t,..

(4.49)

(4.50)

Method

of

lines 183

By substituting

klU~

and

k2U2

from the Riccati relations (4.49) into the

discretized boundary condition in

the

last of (4.47), we see that the

position of

the

interface at t =

t,.

is

the

root x =

Sn

of the equation

,

~l(X,

t,.)-w(x)

{R(X)~2(X,

t,.)+z(x)}

U(x)

+ A (x,

t,.)(X-;;-1)-~3(X'

t,.)=O. (4.51)

Once s(t,.) is known, the functions Uj are obtained by backward integra-

tion of

dU2

1

dx

=k

2

(

RU

2+

Z

),

(4.52)

from the second of (4.49), and

dcfJl

=

Cl

UcfJl

+W

-VI

+ F

dx St I>

with

A..

(

)=~l(S)-W(S)

'P'l

S U(S) ,

b

1

<X<S,

from the first equations of (4.47) and (4.49), and Ul(X) follows from

the

first of (4.49).

Meyer

(1978c) suggested integration of the differential equations by

the simple trapezoidal rule coupled with linear interpolation between

nodal values

to

locate the x =

s,.

of (4.51).

He

discussed the relative

merits of this and more elaborate methods.

Meyer

(1977a,b) treated two- and three-dimensional problems by

coupling his one-dimensional algorithms with the locally one-dimensional

methods familiar in fixed boundary problems, e.g. fractional-step methods

in alternating directions (see, for example, Williams 1979).

Thus, Meyer

(1977a) considered the problem defined by

au

v .

{K(x,

t)Vu}+

a(x,

t) .

Vu

+b(x,

t)u

-c(x,

t)ai=

f(x, t) (4.53)

for x:::

(XI>

X2)

in the domain

D(t),

t>

to,

with conditions u =

a(x,

t), for

points

x on those parts of the boundary aD

l

(t),

t>

to,

which coincide with

the

coordinate axes, and u = go(x, t) for points x on the free boundary

aD

2

(t), together with the second free boundary condition

(4.54)

184 Front-tracking methods

The initial value u(x,

0)

= uo(x) within

the

initial domain D(to)

is

also

given.

One

half of the equation (4.53)

is

discretized in time and integrated

in alternate directions, i.e. first

LlU=.!...(K(X,

t)

aU)+al(X'

t)~+!b(x,

t)u

aXl aXl aXl

au

-!c(x,

t)

at

=!/(x,

t),

is

integrated from

t"

to

t,,+!

with (4.54) replaced by au/ax

l

=

gl(X, t) - Al

asl/at.

The initial value

is

u(x, t,,)

on

D(t,,),

X2

is

considered

as a free parameter and the other boundary conditions are satisfied. This

one-dimensional problem

is

solved by the method of lines and invariant

imbedding (Meyer

1977a)

to

yield values of u(x,

t,,+!)

and D(t,,+!) which

then serve

as

initial values for

the

second, one-dimensional problem

defined by

a (

au)

au

~u = -

K(x,

t) - + a2(x, t)

-+!b(x,

t)u

aX2

au

-!c(x,

t) at

=!/(x,

t),

to

be discretized in time and integrated over

the

interval

t,,+~

to

t,,+h

and

where

Xl

is

now a free parameter. Equation (4.55)

is

solved N times for N

different values

of

the parameter

X2

=

x~,

k = 1, 2,

...

,N.

These are

uncoupled equations and can

be

solved simultaneously, followed by

the

solution of N equations (4.56) for N different values of

Xl

=

x~,

k =

1,2,

...

, N.

Meyer (1977 a,b) pointed

out

that

the

disadvantage of this alternating

direction method

is

its inability

to

track free boundaries which move

essentially along one of the coordinate axes, so that

the

boundary cannot

be

described as an invertible function with respect

to

two orthogonal

axes. Point sources may also lead

to

difficulties.

Meyer

(1978a)

therefore proposed a second method based

on

the

method of lines combined with an

SOR

type iteration, where

the

sub-

sidiary, one-dimensional problems are again solved by invariant imbed-

ding.

The

free boundary

at

any given time

is

visualized as a curve

or

surface over a base domain, e.g. Z =

z(x,

y,

t), r = r(

(J,

C/>,

t) in three

dimensions

or

y = y(x, t), r =

«(J,

t) in two dimensions.

The

spacial vari-

ables of the base domain as well as time are then discretized and all

derivatives with respect

to

these variables are replaced by appropriate

difference ratios. A system of coupled free-boundary problems for

second-order equations in

the

remaining space variable results. This

system

is

solved by cycling through it using a Gauss-Seidel

or

SOR

Method

of

lines

185

algorithm and using invariant imbedding to solve the one-dimensional

free-boundary problem on each line.

As

an example, consider the problem defined by the following system

au

V

2

u

--

= f(x,

y,

t)

at

in the time dependent

domain

D(t),

(4.57)

u =

a(x,

y,

t)

on

the fixed boundary

aD

1

(t), (4.58)

g1(

x,

y,

t,

u,

::

'

::)

=

o}

au

as

(

)

on the free boundary

aDit).

(4.59)

g2

x,

y,

t,

u,;-,-

=0

on at

It

is assumed that

the

free boundary

aD

2

(t) can

be

expressed

as

y = s(x, t)

and

that

any line x = constant cuts

aDit)

at

most once.

In

contrast

to

the

fractional step method, however, it

is

not

necessary for

aDit)

to

be

expressible as

the

inverse function x =

u(y,

t), say. After discretizing time

and

x

the

following equations are to

be

solved

at

time level

t..:

~~

1 1

oy2 + (ax)2

(~+1

-

2~

+

~-1)

-

at

(~-

~.n-1)

=

f(~,

y,

t..)

in D(t..),

~

(y) =

a(~,

y,

t..)

on

aD

1

(t..),

j=1,2

(4.60)

where

a~

=

{(~+1

-

~-1)(S;+1

-

S;-1)_

a~

(S;)}/{1

+

(S;+1

-

S;_1)2}!.

an

2ax

ay

2ax

(4.61)

This

is

a multi-point system of ordinary differential equations with free

boundary points.

A'

successive over-relaxation scheme of solution, start-

ing with an initial guess

(u?,

s?)f:.1>

proceeds by solving in the

kth

iteration

the

equations (4.60) for ii; where

~+1

=

u~;l,

~-1

=

U~-1>

Si+1

=

s~;l,

S;-1

=

S~-1>

i.e. new values are used as soon as they are available in

the

usual

SOR

fashion. Invariant embedding can conveniently be used to

solve

the

sequence of one-dimensional problems and the location of the

free boundary

is

essentially uncoupled from the solution of the second-

order equation in (4.60). Once

(ii;,

sn

are found

the

solution for the next

line

is

obtained by incorporating

the

values

u~

=

U~-1

+ w(ii; -

U~-1),

where w

is

a relaxation factor chosen by experience to maximize the

186

Front-tracking methods

speed of convergence with hope of an improvement from one time level

to

the

next.

The

attraction of using

the

invariant imbedding method of solution, as

in the one-dimensional problems discussed earlier in this section,

is

that it

is

little affected by

the

complicated form of the free-boundary condition

involving (4.61).

Meyer

(1977b) dealt similarly with an elliptic problem in radial coordi-

nates.

He

mentions matters of the existence and uniqueness of his

solutions and anticipates satisfactory convergence

on

the evidence of

computed examples.

Meyer

(1978a) solved a Stefan and an ablation problem by

the

method

of lines with SOR. Furzeland

(1979b) produced a general-purpose com-

puter package, incorporating Meyer's ideas, for one-dimensional, non-

linear moving-boundary problems and

t~ted

it

on a large variety of

problems. George and Damle (1975) considered non-uniform initial

temperatures.

The

method of lines with a Riccati transformation was applied

to

the

binary alloy solidification problem by Meyer (1981).

5.

Front-fixing methods

5.1. One-dbnensional problems

.AN

alternative

to

tracking

the

moving front

is

to

fix

it by a suitable choice

of

new space coordinates.

For

example,

in

the

simple one-dimensional

melting problem specified by equations (1.21-6),

the

transformation

~

= x/s(t)

(S.l)

fixes

the

melting boundary

at

~

= 1 for all

t.

By using

the

standard

relationships

au

1

au

ax

= s(t)

af

au)

au

a~

(au) x

ds

au

(au)

at

x =

a~

at

+

at

f = - {s(tW

dt

a~

+

at

f

to

transform from u(x, t)

to

u(~,

t),

the

heat

equation (1.21) can

be

written

0<~<1,

where s = s(t),

and

(1.26) becomes, for A = 1,

lau

ds

-;

a~=dt'

~=1,

t>O.

t>O,

(S.2)

(S.3)

The

transformation

(S.l)

was proposed by

Landau

(19S0) and first

applied'

to

a finite-difference scheme by

Crank

(19S7a). Lotkin (1960)

used unequal intervals

in

~

and t and divided differences

to

obtain

an

economic improvement in accuracy. Temperature-dependent thermal

properties were included by Citron (1962) and Mastanaiah (1976). Ferriss

and Hill (1974) solved

the

implicit oxygen diffusion with absorption

problem (see §1.3.10) by using

the

same transformation. Some

of

their

results

are

included

in

Table 6.3 for convenience. Dahmardah and

Mayers (1983) used

the

variable x/s(t)

in

a Fourier-series solution of

the

oxygen problem which supports

the

results

of

Hansen

and

Hougaard

(1974)

in

Tables 4.1

and

4.2 except in

the

final stages. Further examples

188

Front-fixing methods

are discussed in Hoffmann (1977, Volume III, pp. 4, 49, 91).

Hohn

(1978) proved second-order convergence of an explicit finite-difference

scheme based

on

(5.1).

If

a fully implicit numerical scheme is used

to

solve (5.2) together with

the

necessary conditions

on

a fixed boundary,

~

= 0, and say u = 0

on

~

= 1, the additional condition (5.3)

on

~

= 1 can

be

incorporated in

the

way eqn (5.15) is treated later in a two-phase problem. Alternatively,

and

necessarily for

an

implicit type of moving boundary condition, in which

ds/dt does

not

appear, some estimated value of s must

be

iteratively

improved in each time step. Thus

the

oxygen diffusion problem has

conditions

u = 0,

au/ax

= 0

on

the

moving boundary and so Ferriss and Hill

(1974), having introduced a variable

cf>

=

u+t,

used

acf>/a~=o

on

~=

1 in

order

to

solve

the

boundary-value problem equivalent to (5.2) for two

estimated values of

s,

denoted

by

SI

and

S2.

Since

cf>

= t

on

~

= 1, s must

satisfy

F(s)

=(</>/t)-l

=0

on

~

= 1. Thus

if

F(SI)

and

F(S2)

are

evaluated

from

the

solution of

the

equivalent of (5.2)

the

method of Regula-Falsi

yields an improved estimate

S3

given

by

S3

= {SIF(SI) - S2F(S2)}/{F(S2) - F(SI)},

and

so on. When advancing

the

time

to

the

next step,

the

initial estimate

for

s is taken

to

be

the

final value

at

the

previous time level.

In attempting

to

obtain quadratic convergence, Baumeister et

al.

(1980)

developed a Newton-type iterative scheme. Following

the

ideas of

Bon-

nerot and Jamet (1974)

and

Fasano

and

Primicerio (1977),

the

Newton

step requires

the

determination of

the

zero of

an

integral form of

the

moving boundary condition

au/ax

=

-ds/dt,

namely

F(s)=s-b+

r g(T)d-r+

r<t)U(x,

t)dx-

ff(X)dx,

where

s(O)

= b and g and f

are

boundary flux and initial temperature

functions. Although quadratic convergence

is

secured,

the

cost in comput-

ing time is considerable because an auxiliary heat

flow

problem must

be

solved on a given

but

time-variable domain in each step of

the

iteration

in order

to

obtain

the

necessary Frechet derivative.

Nitsche (1980) used

the

boundary-fixing transformation (5.1) coupled

with a new

time

variable,

'T

=

11

S-2(t) dt,

to

obtain a weak formulation of

a single-phase Stefan problem with zero-flux condition

on

the

fixed

boundary. A finite-element discretization was effected

and

used

to

estab-

lish existence

and

regularity results

but

no

numerical solutions

were

obtained. Several generalizations of

the

method, including an application

to

the oxygen diffusion problem, were mentioned together with an

extensive list of references

to

papers

on

finite-element formulations

and

free-boundary problems in general.

One-dimensional problems

189

Furzeland (1980) considered the more general, linear equation with

constant coefficients for the unknown

U;

(x,

t) in each phase i = 1 and 2,

cPu;

aU;

au;

k;-2

+a,-+b;u;-c;-=/;(x,

t),

ax ax at

O<t<t*,

(5.4)

with phase 1 defined as .tl

~x~s(t)

and phase 2 as

s(t)~X~.t2'

where

s(t) denotes the moving boundary. Initial and boundary conditions of the

form

Ul

=

UlO(X),

.tl

~x~s(O)}

U2=U20(X),

s(0)~x~.t2

al

(i,

t)U;

+ a2(i, t)

aU;

= a(i, t),

ax

t = 0,

s(O)

given, (5.5)

O<t<t*,

(5.6)

are allowed for, together with

the

following general form

of

the

moving

boundary conditions:

Ul

=

~

=

f.L(s,

s,

t) }

~

aUl

aU2

au

.)

0

x=s(t),

O<t<t*.

(5.7)

u,

- - -

sst

=

ax' ax'

at'

, ,

The

function F need not

be

a polynomial in its variables, though it

commonly

is, and

the

s term will

be

absent in

an

implicit problem.

In

phase 1, .tl

~x

~s(t),

the

moving boundary can

be

fixed by

the

coordi-

nate

transformation

x-.t

l

~l

=

s(t)-.t

l

'

(5.8)

so

that

x

=.tl

becomes

~l

= 0 and x = s(t)

is

~l

= 1.

The

corresponding

transformation for phase 2

is

(5.9)

Considering phase 1 and dropping the subscript

i = 1 in (5.4)

the

equation

becomes

a

2

u . au

)2

(

)2

au

k72+(s-.tl)(a+c~s)-;+(S-.tl

bu -

s-.t

l

c-

as as at

=

(S-.t

l

)2f(x, t),

0<~<1,

t>O,

(5.10)

after multiplying through by (s - .tl)2 for numerical convenience.

The

usual finite-difference discretization in time

at

t

n

+

1J

= (n + 8) at,

where 8 =

~

for Crank-Nicolson and 8 = 1 for a fully implicit scheme, is

combined with central differences for space derivatives

to

give the

190 Front-fixing methods

tridiagonal system

8[Aj

-

Bj]uj~l+[CS2+8(2Bj

-atS

2

b)]ur+

l

-8[A

j

+

Bj]uj:l

=

(1-8)[B

j

-Aauj-l

+[cS

2

-(1-8)(2B

j

-MS

2

b]uj

+(1-

8)[Aj + B

j

]uj+1-

atS

2

!(Xj, t

n

+1),

(5.11)

where

The

fixing of the moving boundary at

~

= 1 allows

aU;

1

au;

ax =

(s-(;)

a~'

i=1,2,

(5.12)

to

be

expressed as standard finite-difference approximations. Thus,

Furzeland

(1980) examines three of

the

possible choices for approximat-

ing

aUl/a~,

i.e.

(3UN

- 4UN-l +

UN-2)

/2a~,

and

(UN+1

- UN-l)/2

a~,

(5.13)

where

j = N

is

~

= 1 and

UN+l

is

a fictitious point, and compares them

numerically for a test problem of Stefan type.

The

first of (5.13) was

inferior

to

the

other two forms. Moving-boundary conditions with an

explicit relationship for

s such as

u=F

2

(s,

t) }

.

=F

(

aUl

aU2) x

=s(t),

s 1 S,

t,

,

ax ax

~=1,

(5.14)

(5.15)

are easily incorporated into an iterative scheme for

the

solution of (5.11)

by approximating (5.15) at time t

n

+

1J

as

(5.16)

and using

one

of (5.13) for

the

space derivatives. Moving boundary

conditions of implicit type, with

s absent, as in

:~=Fl(i,S,t)

}

or

:=(F2~:~tL2)=0}

(5.17)

F

2

(s,

t,

u1>

U2)

= 0 1

S,

t,

ax'

ax (5.18)

can

be

incorporated by using one of (5.17) as the known boundary

One-dimensional problems

191

condition for (5.11) and estimating s from a regulia falsi solution of the

corresponding eqn (5.18) as Ferriss and Hill (1974) did.

Equation (5.11)

is

solved

to

find u

n

+1.1 with an initial guess

sn+I,O

obtained, for example, by extrapolation from the two previous time

levels. Then an improved estimate

sn+1,1

is

obtained from the moving

boundary conditions, as described above, and the scheme iterated until

convergence

is

achieved

to

a prescribed accuracy (usually

2-3

iterations in

Furzeland's test problems for a precision of 10-

4

).

If

s(t)

= e

l

for a finite

length of time then (5.8) is singular.

The

original eqn (5.1) needs

to

be

solved by using, for example, a Crank-Nicolson scheme on a fixed

domain until

s(t)-e

l

becomes non-zero. This

is

not necessary

at

t=O

as

long as

s(O

at»o in (5.11).

Furzeland

(1979b)

developed a general-purpose computer program

based on

the

above iterative algorithm. Mastaniah (1976) used essentially

the

same algorithm to solve problems with temperature-dependent ther-

mal properties, i.e.

k and c depend on

u.

Thus

k(u)

was evaluated at the

mid points

(j

±t,

n +!) by using Taylor's series to extrapolate u

n

+!

from

previous time levels.

An

alternative way of treating eqn (5.4), suggested in a private com-

munication by C. M. Elliott

to

Furzeland (1979b, 1980),

is

to discretize

only

the

space derivatives and

to

integrate the resulting ordinary differen-

tial equations in time along constant

~

=

~

lines.

For

example, eqn (5.2)

is

approximated by

dUj

1

Uj-I-2Uj

+

Uj+l

S

(Uj+l-

Uj-l)

dt

=

S2

(a~)2

+

~j

s 2

a~

,

j = 0,

1,

...

, N

-1,

0<~<1,

t>O.

(5.19)

Moving boundary conditions containing an explicit relation for S can

be

approximated, for example, by

UN

= 0, s =

-(3UN

-

4UN-1

+

UN-2)/2

a~,

t>O.

(5.20)

An

attraction of this approach

is

that well-established algorithms, e.g. the

Gear

algorithm, for stiff ordinary differential equations such as (5.19) can

be

used. Two readily available implementations are the

NAG

algorithm

D02AJF

or

an algorithm by Aitchison (1975) adapted for sparse mat-

rices.

The

stiff-system algorithm will automatically generate dUj/dt, j =

0,1,

...

, N

-1,

ands,

and

will

produce a solution vector

(uo,

Ul""

UN-h

s)

at required time intervals, provided an explicit boundary condition of type

(5.20) exists. This method

is

readily adaptable to non-linear equations

and moving boundary conditions, and

to

multi-component systems in-

cluding combined heat and mass transfer. Implicit conditions on the

Crank J. Free and Moving Boundary Problems

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