Crank J. Free and Moving Boundary Problems

202

Front-fixing methods

Because of the symmetry of this problem, only the quadrant of the prism

enclosed between

x = 1, Y = 1 and the axes x = Y = 0 need

be

considered,

and there are symmetry conditions

au) au)

-

--

-0

ax

x=o

-

ay

y=o

- .

(5.58)

Crank and Gupta (1975) applied

the

IMM transformation to the Y

variable only in equation (5.54) and obtained

ay

=

_{a

2

u _

aZy

la

y

)-3}

ay

,

at ax

2

au

2

\au

au

(5.59)

using the relationships (5.47) with Y written for

x.

The

IMM form of

the

appropriate Patel expression (see equation (1.56) in §1.3.9) for

the

usual

Stefan boundary condition

au/an = -f3v

n

is

(5.60)

Values of

y are calculated on a (u,

x)

mesh for successive time steps 5t

and

y~j

signifies

y(i

5u, j 5x, k 5t) in

the

usual way. Crank and

Gupta

replaced (5.60) by

the

explicit form

..

YN,i

YNj

_l+.

YN

,i

YN,i-l

k_Uk

'

k+1_

k 1 {

(k

_ k

)2}

5

5t

f3

5x YN,i

YN-l,i

(5.61)

where N

5u = 1.

To

be

consistent with this first-order backward-difference approximation

to

(aylau)-l they approximated (5.59) also explicitly by

k+l

k

(k

k

):'2

k 2 k k

Yij

-

Yij

_

Yij-

Yi-l,i

~

+

Yi-lj-

Yi,j-

Yi+1,i

5t 5u ax

2

(y!<.-

Y'-l

.)2

.

I~J

1

,J

(5.62)

In order to discretize

(aZu/ax

2

)y

at

the

point (i 5u, j 5x) values of u are

required, at three equally spaced values of

x, for which Y =

Yij

= y(i 5u,

j 5x). Crank and

Gupta

(1975) interpolated linearly the values of u

corresponding

to

Yi,i

at

Xj-l

and

Xj+l>

e.g. at

Xj-l

they used

U;+l(Yij-l

- Yi,j) -

U;(Yi+l,i-l

- Yi,i)

u=

.

(5.63)

Yij-l

-

Yi+l,i-l

In order to

take

maximum advantage of

the

symmetry of their problem

they used special interpolation and extrapolation procedures on

the

axes

and on

the

line Y =

x:

They used the one-parameter integral method of

Poots

(1962a) (see §3.5.6)

to

calculate

the

positions of some isotherms

and the interface a short time after

the

start of solidification, and

continued from there with

the

!MM.

Reasonable agreement with earlier

Isotherm migration method

203

TABLE

5.1

Values

of

the y coordinate on the solid-liquid interface

for

fixed values

of

x at various times. Solution starts from the

values taken from the

Poots

(1962a) one-parameter

method at

t=0.0461.

8x=8u=0.1;

8t = 0.0001, {3=

1.561

x

0.0 0.1

0.2

0.3 0.4 0.5 0.6

0.05

0.8125 0.8106

0.8048

0.7940

0.7764 0.7476 0.6904

0.10 0.6979 0.6965

0.6921 0.6836

0.6683

0.6392 0.5606

0.15

0.6157 0.6141 0.6095 0.6000

0.5810 0.5201

0.20

0.5473

0.5453

0.5394 0.5268

0.4789

0.25

0.4865

0.4838 0.4755

0.4567

0.3894

0.30 0.4302 0.4263

0.4146 0.3654

0.35

0.3766

0.3708

0.3534 0.2859

0.40

0.3337

0.3158

0.2623

0.45

0.2816 0.2585 0.1893

0.495

0.2376

0.2056

0.1097

results by Allen and Severn (1962) and Lazaridis (1970) was demon-

strated. Their values for

the

movement of the freezing front are given in

Table 5.1 and their isotherms are shown in Fig. 5.5.

Crank and Gupta (1975) discussed how the fixed boundary conditions

influence the choice of

;PU/

ay

2

or

ifu/ax

2

as the term

to

which the IMM

transformation should best

be

applied. They also indicated how alterna-

tive conditions on the fixed boundaries can

be

handled.

1.0

u-o

u=O.l

0.8

y

t

0.6

u-o

0.4

0.2

1.0

-x

FIG. 5.5. Positions of the isotherms having temperatures u =0.0(0.1)1.0 at t =0.495

204

Front-fixing methods

A similar extension of IMM to two dimensions was used by Durack and

Wendroff (1977).

(ii) Complete transformation

of

variables.

The

interpolation used

to

evaluate

ifu/ax

2

in (5.59) becomes more burdensome in three space

dimensions and, in any case, approximations such as (5.63) are of lower

accuracy than

the

replacement of the other second derivative.

An

alterna-

tive is

to

generalize the IMM for multi-dimensional problems by using a

complete transformation, in which

one

space variable is chosen as

the

new dependent variable and

the

independent variables are the remaining

space coordinates and temperature.

The

complete transformation

is

conveniently performed following, for

example, Boadway's (1976) treatment of fluid flow problems.

The

equa-

tion for heat

flow

in a homogeneous medium in which

the

heat conductiv-

ity

K and specific heat c may be functions of temperature

U,

and p

is

the

density can be written

au

a (

au)

a (

au)

a (

au)

cp-=-

K-

+-

K-

+-

K-

.

at

ax

ax ay

ay

az az

(5.64)

Here, temperature

u is a function of

the

space coordinates (x,

y,

z)

at

a given time i.e.

u =

u(x,

y,

z),

(5.65)

and hence

au

du=-dx+-dy+-dz.

ax ay

az

(5.66)

We introduce two dummy variables

q,

and ,

with

q,

=

q,(x,

y,

z),

,=

'(x,

y,

z)

(5.67)

and

relations equivalent to (5.66) for

dq,

and

d,.

We require new

functions

x = x(q" " u),

y = y(q" " u),

z = z(q" " u)

(5.68)

for which we have

ax

dx=-dq,+-d,+-du

aq,

a,

au'

(5.69)

and

similarly for dy, dz. Then by using (5.69) together with (5.66) and

the

Isotherm migration method

205

corresponding expressions for

dct>

and

dC

we

obtain

ax

~u

au

au)

dx=-

-dx+-dy+-dz

au

x ay az

ax

(act>

act>

act»

+-

-dx+-dy+-dz

act>

ax

ay az

ax

(ac

ac

ac)

+-

-dx+-dy+-dz

,

ac

ax

ay

az

(5.70)

together with two equivalent expressions, one for dy and one for dz.

By collecting

the

dx

terms in the three expressions we obtain

ax ax ax

act>

--

act>

ac

au

ax

ay ay ay

ac

act>

ac

au

ax

(5.71)

az

au

0

-

act>

ac

au

ax

and hence, on solving,

act>

= e

y

az _ az

aY)/A,

ax

ac

au

ac

au

(5.72)

where A

is

the

determinant of

the

matrix in (5.71), and similarly for

a{jax

and

au/ax.

Collecting dy terms yields corresponding expressions for

act>/ay,

a{jay,

au/ay,

and

the

dz

terms give

acf>laz,

a{jaz,

au/az.

Proceeding to second

derivatives we require, for example,

~

(Kau)=.i.

(Kau)

act>

+~

(Kau)

ac+.i.

(Kau)

au

, (5.73)

ax

act>

ax

ac

ax

au

ax ax

into which expressions of type (5.72) are

to

be substituted for all

the

first

derivatives. Similar general expressions can

be

obtained for

a(K

au/ay)/ay

and

a(K

au/az)/az.

Finally, for the IMM we wish to interchange temperature u with one

space variable, say z, i.e.

to

write z =

z(x,

y,

u),

and reference

to

the

last

of (5.68) shows that we must now identify

ct>

with x and C with y in the

above analysis. When this

is done in

the

development of (5.73) and

corresponding expressions for

a(K

au/ay)/ay

and

a(K

au/az)/az

noting

that A

=

az/au

and also that

e:t

=

-(!~)(~i).

(5.74)

206

Front-fixing methods

the

heat eqn (5.64) can

be

written finally as

cp

iJz

=

(iJZ){~

(k

iJZ)

+~

(k

iJZ)}

iJt

iJu

iJx iJx

iJy iJy

_

iJk

_

{iJZ)

~

(k

iJZ)

_

iJz

~

(k

az),

iJu

\ax

iJu

iJx

iJy

iJu

iJy

(5.75)

where for convenience

k denotes

K/(iJz/iJu).

In

the special case

of

K,

c,

p

constant, (5.75) reduces

to

cpiJz

=

[fiJ

2

Z

+

iJ

2

Z)(iJZ)2

+

iJ2

z

{1

+

{iJZ)2

+

{~Z)2}

K

iJt

\ax

2

iJy2

iJu iJu

2

\ax

\ay

-2

iJ::u

~;)~~)

-2

iJ~2:U

~;)(:~)

]~~r2.

(5.76)

An

iterative finite-difference solution of

the

steady-state form

of

this

equation in a three-dimensional free-boundary problem is developed for

eqn (8.70)

in

§8.4.2.

In

a Stefan moving-boundary problem, a three-dimensional extension

of Patel's condition (5.60),

to

be

used in conjunction with (5.76)

and

other

fixed-boundary conditions, is

(5.77)

Turland (1979) derived (5.75) by performing

the

complete transformation

in a manner directly relevant

to

the

present problem.

He

then

regrouped

the

terms in (5.75) by using

the

identities, in subscript notation,

(5.78)

and

(5.79)

where

K/z

u

has

been

reinserted for

k.

When (5.78)

and

(5.79)

are

substituted into (5.75) together with an equivalent grouping of

the

corresponding y terms, (5.75) can

be

written in

the

form

cp

iJz

=

_(iJcI>U

+

iJ<Px

+

iJcI>Y)

, (5.80)

iJt

iJu

iJx

iJy

where

(5.81)

Isotherm migration method

207

Physically, the

4>

terms represent heat

fluxes

and Turland (1979) prefer-

red the form (5.81) because it leads readily to difference schemes in which

the heat content of

the

system

is

conserved, which may not

be

true of

numerical schemes based on (5.75). Conservative difference schemes of

this kind were derived by Potter (1973) for the heat-flow equation in its

usual form, by considering

the

heat balance of a volume element. Turland

(1981) derived (5.80) directly in the same

way

by considering

the

fluxes

into

and

out

of the volume bounded by

the

planes in which u and u + du,

al

and

al

+dab

a2

and a2+da2 are constants where, for generality,

ab

a2,

a3

are taken to

be

curvilinear coordinates and u still denotes

temperature.

The

thermal properties of

the

medium are allowed

to

be

temperature dependent,

but

otherwise

the

medium

is

homogeneous.

In

the

orthogonal curvilinear system the metric is defined by

ds

2

=

hi

dai+

h~

da~+

h;

da;,

and

a volume of coordinate space by

where

(5.82)

(5.83)

(5.84)

The

intention is

to

interchange

the

curvilinear coordinate

al

with temper-

ature

u.

The

functional integral I(u,

u+Au,

a2,

a2+Aa2, a3, a3+Aa3,

t)

is therefore defined by

l

u+AU

[V(U+AU)

1=

pcVdu'=[hVJ::+Au-

hdV'

(u)

(5.85)

where

h =

m(a,)

pc

du'. Thus, correspondingly,

aI

=

JJJpc

az

da2

da

3 du

=-Jpc

au dV'.

at at at

(5.86)

The

partial derivatives of I

and

Z with respect

to

time are taken with u

and u + Au constant, though

the

corresponding value of

al

may change.

The

appropriate heat

flow

equation

is

au

pc-=V·

(KVu)

at

(5.87)

and hence Gauss's theorem coupled with

the

second form of aI/at in

(5.86) yields

aI

#

-=-

KVu

·dS

at '

(5.88)

208

Front-fixing methods

where the integral is taken over

the

surface of

the

volume a v formed by

the intersections of the surfaces

u and u +

au,

a2

and a2 + aa2, a3 and

a3

+ aa3' This volume is

not

a cuboid in shape because its specification

includes the pair of surfaces of constant

u values

rather

than

constant

al'

The

heat flow

per

unit

area

at

a given point, perpendicular

to

planes

of

constant

u,

is

{

1

au A 1 au A 1 au A }

-KVu=-K

--"l+--~+--"3

,

hlaal

h2aa2

h3

aa

3

(5.89)

where

al>

a2,

a3

are unit vectors.

The

element of

area

dS(u)

of

the

constant

u surfaces is

(u)

Vu

dS = hlh2

h

3

aa2

aa3

-1-

,

au aal

where aU/aal is taken

at

constant a2, a3,

and

t.

Thus,

the

net contribution

to

the

heat flux from the sides u and u +

au

of

the

volume element is,

to

first order,

a1(U)

a {

Vu

}

--=

--

(KVu)hlh2h3aa2aa3-1-

au.

at

au au aal

For

the constant

a2

surfaces

()

h h A A A hlh3 A A A

dS

<><,

= 1 3

aal

aa3"2

=

--

aU

aa3"2

au/aal

which, when coupled with (5.89), gives

the

net

contribution

at<><,)

a

{1

au

hlh3}

--=-K-

----

aUaa3aa2'

at

aa2

h2

aa2

aU/aal

(5.90)

(5.91)

The

expression for ato.,)/at follows similarly. Finally,

by

collecting terms

we find

a1

at

u

)

a1(<><')

a1(O<,)

-=--+--+--

at

at at

at

=_[~

(h1h2h3KIVuI2)+K~

(hlh3au/aa2)

au

aU/aa1

aa2

h2

aU/aa1

+ K

~

(hlh2 au/aa

3

)]

au

aa2

aa3'

aa3

h3

aU/aal

(5.92)

To

complete

the

interchange of u and al> which

is

to

become the new

dependent variable, we write

and

Isothenn migration method

209

Then, in the limit

au,

aa2,

aa3

all approaching zero

the

first of (5.86)

becomes

aJ

aZ

-=pc-

aa2aa3au

at at

(5.94)

and

differentiation of (5.84) gives

az/at

= hIh2h3 aal/at. (5.95)

Inserting (5.89) and (5.93), (5.94), and (5.95) in (5.92) we obtain the

partial differential equation for al(a2, a3,

u,

t) which can

be

written in the

conservative form of

(5.80), i.e.

(5.96)

where

<I>(U)

=

h2

h

3

~

{1

+

h!

(aal)\

h!

(aa

l

)2},

hI

aal/aU h2

aa2

h3

aa3

<I>(<<Z)

=

-hIh3

K

aal

<I>(aJ=-

hIh2 K

aa1

.

h2

aa2

h3

aa3

(5.97)

Special cases

of

practical interest can easily

be

extracted from (5.96)

and

(5.97).

For

example,

on

putting

hI

=

h2

=

h3

= 1

and

identifying al

with

Z,

a2 with

X,

and

a3 with y we regain (5.80) and (5.81) which were

seen

to

be

equivalent

to

(5.75).

The

usual heat-flow equation in

one

space

variable can

be

written in

the

form

pc au

=.!.i.

(Krp

au),

at

r

P

ar ar

where for plane geometry p = 0 and r =

X,

and for cylindrical

and

spheri-

cal symmetry r is

the

radial coordinate and p = 1

or

2 respectively.

For

plane geometry hI =

h2

= h3' for cylindrical coordinates hI = 1,

h2

=

r,

h3

= 1, and for spherical polar coordinates

hI

= 1,

h2

=

r,

h3

= r sin

6.

Thus

IMM

expressions in three dimensions for the plane, cylindrical, and

spherical geometries are easily deduced from

(5.96) and (5.97) by ap-

propriate choice of

hi>

h2'

h

3

.

In

particular, the IMM expressions in

one

space variable become

ar

1 a (Kr

P

)

pc

at

= - r

P

au

ar/au '

with p

= 0, 1, 2 as above.

(5.98)

The

two-dimensional case in spherical polar coordinates (of interest for

example in studies of

the

melt-pool growth in a reactor substrate (Turland

210

Front-fixing methods

1979)) is the case of

h1

= 1,

h2

=

r,

h3

= r sin

6,

and

in (5.96), (5.97),

a1

=

r,

a2

=

6,

and

a3

terms are ignored. Thus (5.96) becomes

-per

2

sin 6

ar

=~

[Kr

2

Sin

6 {1 +.!.lar)2}]

_~

(K

sin 6

ar).

(5.99)

at

au

ar/au

r2

\a6

a6

In

order to illustrate how equation (5.96), written in terms of the

heat

fluxes

<1>,

leads to a conservative form

of

difference equation, we consider

the

discretization of the one-dimensional plane case of (5.96) and its

derivation as presented above. Thus

the

one-dimensional analogue

of

(5.85) may

be

written (Turland 1981)

I(ul>

U2,

t) =

1:

2

pcx

dUo

and hence, using (5.98) with p = 0,

aI f U

2

ax

[U

2

a ( K )

[K

]U

2

at

= J

Ul

pea;

du

= - ,au

\ax/au

du

= -

ax/au

u,'

(5.100)

Defining

the

heat

content

per

unit area,

H,

by

H(xl>

~,

t) =

fx~

h

dx

and

h =

JO(x)

pe du,

then

the

equivalent of (5.85) is

(5.101)

where

Xl>

h1

correspond to

Ul>

etc.

It

follows from (5.101) that

if

Xl>

U1

and

X2,

U2 are fixed (being

the

sides of a volume element) then conserva-

tion of

the

quantity I implies conservation of total heat.

In order now

to

set up difference equations for numerical solution we

can consider

the

set

of isotherms llt =

Uo

+ i

au

for

o:s;;

i:s;;

j say. Then use

of (5.100) gives

and

aI [ K

]U;-t

at

(U!,

llj-!,

t)

= -

ax/au

U!·

We approximate by using

I(llt-!> llt+!,

t)

=

peiXj

au

=

n}

-K)

_ -Ki+t

au

*

--

- lLi+!,

ax/au

i+!

Xj+1

-

Xj

and therefore (5.102) and (5.104) become

an

=

14'+!

_

IL1'-J,.

dXj

=..!..

{Ki-!

Ki+!}

at

dt

PC;

Xj

-

Xj-l

Xj+l

-

Xj

(5.102)

(5.103)

(5.104)

(5.105)

Isotherm migration method

211

By summing the first

of

(5.105) for

l~i~j-l,

the difference

fonn

of

(5.103) is

d

j-1

~

r,-

* *

dt

.'-'

i - /LH - /L!,

.=1

(5.106)

which is satisfied exactly

and

so

the

method is conservative.

Turland (1981) referred

to

numerical experiments based

on

a simple,

explicit replacement of

dx,/dt in

the

second of (5.105). Use of an implicit

scheme, however, has

the

usual advantage of removing

the

stability

restriction on

the

time step.

He

proposed the linearized scheme in the

time interval

t

n

to t

n

+

1

n+1

n

Xi

-Xi

_i{

n n

+bn

n+

n n + * n+1+b* n+1+ *

n+1}

+1

n

-2

ai

x

l-1

IXi

Ci

X

i+1 a

i

x

i-1

iXi

Ci

X

i+1,

t

n

-t

where

a;=K;-!/A"

G=K;+t/A;,

bi=-(a;+G),

Ai

= pc,

(x,

-

X,-1)(X,+1-

Xi),

(5.107)

(5.108)

(5.109)

(5.110)

The

values of

x1

calculated from the predictor fonnula (5.110) are

inserted into (5.109) and (5.108) to evaluate the coefficients

a1,

bt,

ct

in

(5.107). With appropriate discretization of boundary conditions the val-

ues of

Xf+1 are then

the

solution of a linear, tridiagonal system of

equations of type (5.107).

The

method is stable for all time steps and

the

criterion

Ixf+1-

x!1 < E for all i can

be

used

to

choose

the

time step.

Turland (1982) solved a number of test problems satisfactorily.

(iii) Isotherm migration along orthogonal

flow

lines. Crank and Crowley

(1978) described a novel approach

to

the

solution of transient heat-flow

problems in two dimensions.

The

movements of isotherms along

or-

thogonal flow lines are tracked in successive small intervals of time by

solving a locally one-dimensional IMM

fonn

of a radial heat equation.

The

changing shape and orientation of the orthogonal system are catered

for by a geometric procedure.

The

procedure is analogous to

the

Lagran-

gian fonnulation of fluid flow, in which the motion of particular particles

of

fluid is calculated,

rather

than

the

velocity distribution at fixed points in

space.

In

non-dimensional fonn, the equation of heat flow in cylindrical polar

coordinates

(r,6) may

be

written

au a

2

u 1 au 1 a

2

u

-=-+--+--

at

ar2

r

ar

r2

a6

2

,

(5.111)

Crank J. Free and Moving Boundary Problems

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