Crank J. Free and Moving Boundary Problems

242

Fixed-domain methods

Crowley and Ockendon (1977) solved a similar problem with body

heating in modelling a thermal switch. They obtained an asymptotic

solution of an integral equation for large times and enthalpy-based

numerical results.

6.2.7. Other conseroation forms and weak solutions

The

essential features of

the

enthalpy method, introduced in §6.2, are

that the problem

is

reformulated in terms of a single differential equation

which

is

solved over a fixed domain

to

obtain values of temperature and

the enthalpy function.

The

position of

the

moving boundary emerges as a

by-product of

the

computation.

In

heat-How problems, enthalpy has a

real physical significance and the reformulated problem

is

in a heat

conservation form. This means

that

a boundary jump condition such as

(6.4)

is

automatically satisfied. Even in problems where the concept of

enthalpy has no direct relevance it

is

sometimes possible

to

secure

the

advantages of computing in a fixed domain by rewriting

the

equations in a

form in which some essential physical quantity is conserved and which

permits a weak solution.

The

following are some examples.

(i)

Cases

of

zero 'specific heat'.

In

the

classical Stefan problem, specific

heat

c and conductivity K are both positive and

the

enthalpy H

is

a

strictly monotone function of temperature, while temperature is a single-

valued function of

H. Crowley (1979) formulated certain problems by

introducing a 'generalized enthalpy' function incorporating a zero 'specific

heat' in one

or

both phases.

An

example of the latter case

is

the

electro-chemical machining problem described by McGeough and Ras-

mussen (1974) (see §2.12.3). Their quasi-steady model

is defined by the

equation

~<f>=o,

(6.66)

where

<f>

is

the electric potential, with the boundary conditions at the

anode

<f>

= 0,

ds

M-=-V<f>·

Vs

dt

and

on

the cathode, r =

a,

on

<f>=-V<o.

s(r, t) = 0,

M>O,

(6.67)

(6.68)

This

is

a one-phase problem governed by an elliptic equation and a Stefan

boundary condition with non-zero latent heat, M.

In

order to formulate

the equations

on

a fixed domain, the anode can

be

regarded as a region

where

<f>

=0,

H=M,

and solutions are required of

(6.69)

Enthalpy method

243

with boundary conditions

c(>

=

-v,

H = 0

on

r =

a,

and initially

cf>=

-v,

H=O,

t=O,

outside s(r, t) = 0,

and

cf>

= 0,

H=M,

t=O,

within s(r, t) =

O.

The

conservative form of (6.69) across a surface of discontinuity s(r, t) =

0, obtained by integrating over a small volume,

is

[H]ds/dt=M

ds/dt =

[Vcf>

• Vs],

(6.70)

where [ ] denotes the jump between the sides where

H~O

and

H;;:.:M

respectively. Since

V2c(>

= 0 inside

the

anode

it

follows

that

(6.67)

is

satisfied by a solution of (6.69). Crowley (1979) obtained an implicit,

finite-difference solution for this problem by expressing (6.69) in cylindri-

cal polar coordinates.

Her

results are compared with those of Christian-

sen and Rasmussen (1976) and Elliott (1980) in Table 8.25.

The

flow of an incompressible fluid in a Hele-Shaw cell (Richardson

1972) (see §2.12.2) can be formulated in the same way (see Figs. 2.19,

20).

In

Crowley'S second problem the 'specific heat'

is

zero only in one

phase in which, therefore, an elliptic differential eqution holds, but we

have a parabolic equation in

the

other

phase. This is the situation for

saturated/unsaturated

flow

in a porous medium (Hornung 1978).

In

one

region the flow

is

unsaturated,

that

is

the medium is only partially filled

with fluid and its fluid content may change.

In

a second region the

medium is saturated and no more fluid can

be

added.

In

terms of a

velocity potential,

C(>,

the

relevant equation

is

(6.71)

where H measures

the

air content of

the

medium. Hornung (1978) adopts

the

'generalized-enthalpy relationships',

c(>

>0,

c(><0,

(6.72)

and

cf>

can

be

regarded as a generalized temperature. A sample problem

is

solved by Crowley (1979).

In

the

same paper, the uniqueness

of

the

weak solutions for these cases

in which

the

'specific heat' may vanish was proved by necessary exten-

sions

of

Oleinik's (1960) methods for the classical Stefan problem.

It

is

because zero specific heat means that temperature

is

no

longer a singie-

valued function of enthalpy

that

Oleinik's proof needs modification.

Crowley (1979) introduced a function strictly monotone in both

H and

C(>,

defined by

F(

cf>,

H)

= H +

'Ycf>

where

'Y

is

any positive constant, into the

244

Fixed-domain methods

uniqueness proof. This has the additional consequence that the finite-

difference scheme used to solve (6.69), for example, cannot be explicit

but

must be implicit. Thus, an implicit scheme of the form

H(x,

t + 8t) -

H(x,

t) =

8t{1-'-

V

2

cf>(x,

t + 8t) +

(I-I-'-)

V

2

cf>(x,

t)}

for ° <

I-'-

.:;;;

1 can evaluate a linear combination of

Hand

cf>

of the form

H

+

"Icf>

used in the uniqueness proof. This follows immediately by writing

the implicit scheme for a one-dimensional problem, for example, as

(H::.+

1

+

21-'-

8tcf>::'+1)

= H::' +

I-'-

8t(cf>::.-t:!-1

+

cf>::'~\)

+

(I-I-'-)

8t(cf>::'-1

+

cf>::'+1

-

2cf>::'),

and solving by successive over-relaxation with

the

right-hand side treated

as known from

the

previous iteration (see solution of (6.51)). From

numerical experiments Crowley (1979) showed

that

the finite-difference

scheme converges but no analytical proof was available because of the

vanishing specific heat. NiezgOdka (1980) and Visintin (1983) estab-

lished the existence of weak solutions for problems with zero specific heat

and examined their stability.

We

note in passing

that

Fasano and

Primicerio

(1979d) proved existence, uniqueness, and continuous depen-

dence theorems for classical solutions for saturated/unsaturated flow in

one space dimension.

Crowley's (1979) approach can

be

adopted in solving the oxygen

diffusion with absorption problem which

is

formulated classically in

equations (1.57-60). By introducing a generalized enthalpy function H

related to oxygen concentration

u by

H={U,

0,

and an absorption term

Q defined by

Q={-I,

0,

the problem can be reformulated as

u>o,

u.:;;;O,

u>o,

u.:;;;O,

O<x<l,

au/ax

=0,

x=O;

u=o,

x=l,

H=u=!(1-x?,

t>o,

t=O.

(6.73)

(6.74)

(6.75)

(6.76)

(6.77)

Integration of (6.75) over a small volume element around the moving

front,

x = s(t), for a time interval 8t

~

0 shows

that

condition (1.59) is

satisfied. Furzeland's (1980) results obtained in this way are included in

Enthalpy method

245

u

FIG. 6.7. Typical equilibrium diagram for a dilute alloy

Table 6.3. They correspond

to

a heat flow problem with zero latent heat

and zero specific

heat

in

the

one

phase where u =

O.

(ii)

Alloy

solidification. Crowleymnd Ockendon (1979) expressed equations

for

the

solidification

of

an alloy in conservation form and obtained numeri-

cal solutions within a fixed domain. They adopted

the

simple linearized

form of equilibrium

or

eutectic diagram (Fig. 6.7) as discussed in §1.3.7.

Denoting by c

the

concentration of impurity and u

the

temperature, then

the

material is in

the

stable liquid

(i

= 1)

or

solid

(i

= 2) phase depending

on

whether c +

k;u

~O

respectively. Assuming

that

the

temperature is

continuous

at

the

phase boundary and that the material

on

either side is

in chemical equilibrium,

there

is

a concentration discontinuity

there

with

Cs

= (kl/k2)CV (6.78)

A one-dimensional solidification process in which the impurity diffuses in

both

phases i = 1, 2 is described by the equations

aG=~(n.

a

G

)

at

ax

'ax'

i = 1, 2,

(6.79)

u

aUt

=

~

(Kt

aUt),

at

ax ax

i =

1,2,

(6.80)

where

D

is

the

mass diffusion coefficient and K the heat conductivity.

The

specific

heat

u of

the

bulk material is taken as unity.

At

the

phase-change interface

the

equilibrium conditions are

Ul

=

U2,

Cl

=

-k1Ul>

C2 =

-k

2

Uz

(6.81)

and

the

conservation

of

heat

and

mass gives

L

dS=_[Kau]2,

[c]ids=_[Dac]2.

dt

ax

1

dt

ax

1

(6.82)

There

are also appropriate conditions

on

the fixed boundaries.

An

246

Fixed-domain methods

enthalpy form which matches the first of (6.82) is

H-{U,

u+L,

u+c/k

1

<0,

u+c/k

2

>0.

(6.83)

Making the usual assumption for dilute solutions,

that

the concentration

of impurity is proportional

to

its chemical potential

or

activity

v,

we have

i=1,2,

(6.84)

where

b is some constant.

The

interface condition (6.78) becomes

Vi

=

V2

and so V, as well as u,

is

continuous across

the

phase boundary u + bv =

O.

Since

~

= constant in this formulation, V can

be

defined everywhere by

(6.84), though in a more general case this relation only holds

on

the phase

boundary.

Thus, the chemical activity and concentration in

the

mass transfer

problem are

the

direct analogues

of

temperature and enthalpy in the

heat

transfer and (6.79) may

be

written in aliform analogous

to

(6.80), i.e.

with

aG

=

~

(D~

aVi)

at

ax

•

ax

'

D;=Dib~,

[v]i=O,

[

c

12 ds =

-[D'

av]2

Jl

dt

ax

l'

i=1,2,

u+bv=O

at x = s(t). Finally, therefore, we can write

the

conservation forms

ac

=~

(D~

av)

at

ax

•

ax

'

aH=~

(K;

au)

at

ax ax

coupled with

C

=bk

1

v,

bk

1

v:S;;c:s;;bk

2

v,

c

=bk

2

v,

H=u,

u+bv<O,

u:s;;H:s;;u+L,

u+bv=O,

H=u+L,

u+bv>O.

(6.85)

(6.86)

Crowley and Ockendon (1979) were unable to prove either existence

or

uniqueness of a weak solution of (6.85) in general. Following Oleinik

(1960) proofs can

be

established for

the

particular case of

Ki

= Ki(u) and

D;=D;(v)

for a

pure

metal where

c=v=O.

Accordingly, Crowley and

Ockendon (1979) proposed a numerical discretization

of

the general case

(6.85), (6.86) with

Ki

and

Di

constant in each phase for simplicity, and

tested their solution by comparison with Rubinstein's (1971) explicit

similarity solutions.

For

the numerical scheme, (6.86) must

be

rewritten

to

give

u,

V in

Enthalpy method

247

terms

of

H,

c,

i.e.

u = H, v =

c/bkt.

H +C/k1

<0,

u=H-L,

v==c/bk

2

,

H+c/k2>L.

(6.87)

The

interval between H

+c/k

1

>0

and H

+C/k2<L

needs special consid-

eration since

u is

not

defined there. Since the material cannot exist in

equilibrium in this range, Chalmers (1964) took

the

values

of

H and c to

represent

an

element of material of which a fraction f is solid

and

1 - f is

liquid. This is reminiscent of the procedure adopted by Voller and Cross

(1981a) (see

eqn

(6.55)).

The

vaues of

u,

CS,

CL

for

the

element are

compatible with local equilibrium. Thus we have

H

==

u +

(1-

f)L,

where

Cs

=

-k

1

u,

CL

=

-k

2

u.

Elimination

of

cs,

CL,

f yields a quadratic

equation for

u with

the

solutiQn corresponding to taking

the

positive

square

root

given by

u

==u*

-{k

1

-

H(k

1

- k

2

)/L}+[{k

1

-

H(k

1

-kz)/L}2_4c(k

1

- k

2

)/Ly.

2(k

1

-k

2

)/L

(6.88)

so

that

u = H

on

the

solidus and u = H - L

on

the

liquidus boundary.

The

inverted form

of

(6.86) is therefore (6.87) together with

u=u*,

v

=-u*/b,

H+c/k

2

-L<0<H+c/kt.

(6.89)

where

u*

comes from (6.88).

A simple explicit, finite-difference scheme for (6.85)

if

the

mesh points

(n

-1)

8x, n 8x, (n + 1)

8x

are

in

the

same phase is

H

m+1_Hm+Km

8t

(m

2 m+ m )

n - n n (8X)2 U

n

+1

-

Un

u

n

-1 ,

(6.90)

m+1

m+D,m

8t

(m

2 m+ m )

(691)

Cn = C

n

(8X)2

Vn+1

-

Vn

V

n

-1,

.

where H~ =

H(n

8x, m 8t), m = 0,

...

,M,

n = 1,

...

,N,

etc., and

K~,

D~m

take

the

values appropriate

to

that

phase.

In

order

to

use (6.90) and

(6.91) successfully over the whole

of

a fixed domain including

both

phases, we require

them

to

conserve

heat

and mass across

the

phase

boundary. Integration

of

(6.90) shows

that

the

change in heat content in

an

interval

8t

is given by

N-1 N-1

8x

L

(H~+l-

H~)

= (8t/8x) L

{K~(u~+1

-

u~)-

K~(u~-

U~-l)}

1 1

N-1

= (8t/8x) L

(K~-

~+1)(U~+1

-

U~)

+

K;:':(u;:':-

U;:':-l)

-

Kl'(ul'-

u;;').

1

248

Fixed-domain methods

Thus, there will

be

a heat loss

or

gain

at

the

phase boundary unless

K;:'

=

K;:'+l

for n = 1,

...

, N

-1,

which implies

Kl

= K

2

•

Modified forms

of

K and

D'

are therefore introduced

near

the

phase boundary, the

position

of

which corresponds

to

u +

bv

= 0

and

is found by linear

extrapolation

on

the

values of u +

bv

at

the

nearby mesh points.

If

s(mcSt)=ncSx+hh

0<h1<cSx, with

the

solid phase in

x<s(t),

the ex-

pressions for the space derivatives of

u in (6.85) are replaced by

where

K*

is

the

modified heat conductivity. Similar expressions are used

for

the

space derivatives

of

v.

Crowley and Ockendon (1979) secured

the

best agreement with

the

appropriate analytical solution by using

K*

=

Kl

h1/cSx

+

Kz(l

-

h1/cSx)

and

similarly for D'*.

At

a mesh point where

u;:'+

bv;:'

= 0,

the

equation

(6.94)

is used.

The

final, finite-difference scheme provides a conservative al-

gorithm

on

a fixed domain which avoids explicit application of

the

classical conditions (6.82) at the phase boundary.

As

a test example, Crowley

and

Ockendon (1979) obtained numerical

solutions for

the

solidification of a block

of

molten alloy, with uniform

temperature and concentration distributions,

Uo=

1,

co=O.l,

occupying

the

region

x>

O.

At

time t = 0,

the

temperature at

the

face x = 0 is

lowered

to

u =

-1

and they assumed

there

is no mass flux

out

of

the

material. They obtained satisfactory agreement, for chosen values of

the

parameters, with Rubinstein's (1971) analytical solution, which is

Cl

=

-k1Us,

erfc{x/2(D2t)~}

c2=0.1-(0.1+k

2

u

s)

erfc{{3ID1} ,

1

(1

) erf{x/2(Klt)!}

Ul

= - + +

Us

erf{{3/Ki}

,

- 1 (

-1)

erfc{x/2(K2t)!}

U2

- + U

3

erfc{f3/

Iq} ,

s(t) =

2{3t!,

(6.95)

Enthalpy method

249

where

(3

is

the

root of

L(37T!

erf«(3/Kl)erfc«(3/Iq)exp«(32/K

1

+

(32/K~

-

K~rfc«(3/Iq)exp«(32/K2)

+ Iqerfc«(3/KVexp«(32/K1)

= {Ki

erfc«(3/Iq)exp«(32/K~

+

K~

erfc«(3/Kf)exp«(32/K

1

)}us, (6.96)

and where

Us,

the

constant phase-change temperature, is given by

O.lD!

(6.97)

The

error

in

the

position

of

the

phase boundary obtained from the

finite-difference scheme compared with

the

analytical solution (6.95)

is

less

than 5

per

cent. Slight oscillations in

the

numerical solution are less

noticeable

the

greater

the

number

of

mesh points.

A more general enthalpy-type formulation by Wilson

et

al.

(1982)

removes

the

restriction that the heat capacities in the liquid

and

solid

phases must

be

equal

and

allows more general boundary conditions

to

be

considered. They give a number of numerical finite-difference solutions

and some comparisons with Rubinstein's solution (1971).

The

possible

appearance of a mushy region just in front of the interface is mentioned.

(iii) Flame-front problems. Sometimes

the

equations in

one

part

of

a fixed

domain are special cases

of

the

equations in the remainder of

the

domain.

In these cases,

the

latter equations are valid for

the

whole domain

but

when written in conservation form they may not satisfy

the

jump

condi-

tion

on

the

moving boundary, which forms

part

of the classical formula-

tion. Crowley (1981) handles this situation by adding a term

to

the

differential equation which vanishes everywhere except

on

the moving

boundary where it has

the

required jump.

The

context is an enthalpy

formulation of Buckmaster's (1979) model of flame fronts, in particular a

steady-state flame tip for laminar, pre-mixed flames with large activation

energy.

The

reaction region is considered

to

be

confined

to

a flame sheet

y = h(x), which divides two regions where the simple diffusion

and

heat

conduction equations are satisfied. Jump conditions hold

at

the

flame

sheet.

The

governing equations considered are

(6.98)

(6.99)

T=l+T"",

y

>h(x).

(6.100)

250

Fixed-domain methods

The

conditions

on

the flame sheet, y = h(x), are

T=1+Too,

[<1>]

= 0 =

[T],

[

a<l>]

[all

{<I>Y=h}

ay

=-A

ayJ=Aex

p

2(1+T.Y

,

(6.101)

(6.102)

(6.103)

where [ ] denotes

the

jump from

one

side of

the

front

to

the

other

as

usual.

The

other

boundary conditions are

aTlay = 0 =

a<l>lay,

y=O,

(6.104)

and

on

x

=0

we have

T(O,

y) = T .. +exp(y - h

o

),}

<I>

= - A (y - ho)exp(y - h

o

),

O<y<h

o

,

x=O

(6.105)

T=1+T

.. ,

<I>

= 0,

y>ho,

x=O.

(6.106)

In

these expressions T denotes temperature,

<I>

measures the sum

of

temperature

and

concentration

of

fuel related

to

its value far behind the

flame, and A

=

(1-1/L)6,

where 6

is

the activation energy and L

the

Lewis number.

It

is

obvious

on

inspection

that

equations (6.100) are simply the

constant-temperature forms of (6.98)

and

(6.99), which are therefore

valid over the whole domain except perhaps

on

y = h(x).

The

conserva-

tion form of (6.99)

on

the

flame front can

be

obtained by integrating

(6.99) with respect to y across the front, having replaced

a<l>lax

by

(a<l>lay)(dhldx).

It

is

[<1>]

dhldx =

-[a<l>lay

+ A aTlay], which shows that

the

first of the jump conditions in each of (6.102) and (6.103) are satisfied.

On

the other hand,

the

conservative form

of

(6.98) gives [T] dhldx =

[aT lay],

which does not agree with the second of each

of

(6.102) and

(6.103). Crowley (1981) therefore modifies equation (6.98) by the addi-

tion of a term which vanishes everywhere except

on

y = h(x), where it has

the jump required in

the

second

of

(6.103). When (6.98)

is

rewritten as

aT

a {H(1 ) G

<l>y=h

)}

a

2

T

-+-

+T

..

-T

exp

(1

T.)2

=-2'

ax

ay

+..

ay

where

H(u)

= 0, u

~O

and

H(u)

=

1,

u

>0,

the

conservation form at the

discontinuity gives

[T]:~

-

[H(1+T

.. -

T)eXPG(:~~

..

f)]

=-[:~,

which at T = 1 + T .. implies the desired jump conditions in (6.102) and

(6.103).

Enthalpy method

251

For

numerical Sblution, the conditions T = 1 +

Too,

e(>

==

0 on y

==.e

for

some sufficiently large

.e>

ho

are added and then the governing equations

in conservation form, to

be

solved over the whole of the fixed domain,

O~y

~.e,

are

aT

a { (

C(>Y=h

)}

a

2

T

ax

+ ay

H(l

+ Too-T)eXPU(l +

Too?

=

ay2'

ae(>

a

2

T a

2

e(>

-=>..-+-

ax ay2

ay2'

together with

the

conditions (6.104-106).

(6.107)

(6.108)

A convenient numerical procedure, given

e(>,

T at all mesh points j

8y,

i8x,

O<y<.e,

is

(i)

find

h(i

8x) where T = 1 +

Too

by quadratic extrapolation on T from

the region

y<h(i8x),

i.e.

T<l+Too;

(ii) evaluate

e(>y=h

==

e(>~

by quadratic extrapolation from both sides and

taking

the

mean;

(iii) solve (6.107) explicitly using

T}<l+Too,

(iv)

solve (6.108) as

-I..~+1

=

-I..~+~

{(-I..!+1_

2-1..~+1-

-I..~+1)

'l"J 'l"J

2(8y)2

'l"J+1

'l"J

'l"J-1

+

>"(T}!t

-2T;+1+

T}~D+(e(>}+1

-2e(>}+e(>}-1)

+>"(T}+1-

2T

}+T}-1)},

using successive over-relaxation.

Typical profiles of the flame sheet are shown in Fig. 6.8. Crowley

(1981) was not able

to

quote existence and uniqueness proofs of weak

solutions for the coupled systems of equations nor

to

prove the con-

vergence of the numerical scheme just outlined.

In

another paper Crowley

(1980) suggests enthalpy treatments of flame fronts in solid fuel and a

quite different problem relating

to

the growth of a crystal as it

is

slowly

withdrawn from a bath of its own melt. The movement of an air-liquid

interface has

to

be

considered as well as that of the solidification front,

with the additional complication that a corner exists all the time in the

interface with air.

Crank J. Free and Moving Boundary Problems

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